us/s// 


MIT/WHOI  97-10 


Massachusetts  Institute  of  Technology 
Woods  Hole  Oceanographic  Institution 


Joint  Program 
in  Oceanography/ 
Applied  Ocean  Science 
and  Engineering 


1930 


DOCTORAL  DISSERTATION 

Dynamics  of  Nonlinear  Cross-Eguatorial  Flow 

in  the  Deep  Ocean 


MIT/WHOI 

97-10 


Dynamics  of  Nonlinear  Cross-Equatorial  Flow 
in  the  Deep  Ocean 


by 

Christopher  A.  Edwards 


Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02139 

and 

Woods  Hole  Oceanographic  Institution 
Woods  Hole,  Massachusetts  02543 

December  1996 

DOCTORAL  DISSERTATION 

Funding  was  provided  by  the  National  Science  Foundation  and  an  Office  of  Naval  Research 
ASSERT  grant  through  the  Massachusetts  Institute  of  Technology.  Computational  resources  were 
also  received  through  the  SCOUT  project  at  the  Massachusetts  Institute  of  Technology. 


Reproduction  in  whole  or  in  part  is  permitted  for  any  purpose  of  the  United  States  Government, 
This  thesis  should  be  cited  as:  Christopher  A.  Edwards,  1996,  D3niamics  of  Nonlinear 
Cross-Equatorial  Flow  in  the  Deep  Ocean.  Ph.D.  Thesis.  MIT/WHOI,  97-10. 


Approved  for  publication;  distribution  unlimited. 


Approved  for  Distribution: 


Philip  L.  Richardson,  Chair 
Department  of  Physical  Oceanography 


V\03\C.vb^  'L.  KAtrCvjCSXT 

Marcia  K.  McNutt 
MIT  Director  of  Joint  Program 


IffIC 


1x4 


D  r> 


WHOI  Dean  of  Graduate 
Studies 


DYNAMICS  OF  NONLINEAR  CROSS-EQUATORIAL 
FLOW  IN  THE  DEEP  OCEAN 

by 

Christopher  A.  Edwards 

B.S.  Haverford  College 
(1988) 

Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

Doctor  of  Philosophy 

at  the 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 

and  the 

WOODS  HOLE  OCEANOGRAPHIC  INSTITUTION 
December  1996 


Professor 


Massachusetts  Institute  of  Technology 
Thesis  Supervisor 

Accepted  by .  . . . . 

Paola  Malanotte-Rizzoli 
Chairman,  Joint  Committee  for  Physical  Oceanography 
Massachusetts  Institute  of  Technology 
Woods  Hole  Oceanographic  Institution 


DYNAMICS  OF  NONLINEAR  CROSS-EQUATORIAL  FLOW  IN  THE 

DEEP  OCEAN 

by 

Christopher  A.  Edwards 


Submitted  in  partial  fulfillment  of  the  requirements  for  the  degree  of 
Doctor  of  Philosophy  at  the  Massachusetts  Institute  of  Technology 
and  the  Woods  Hole  Oceanographic  Institution 
December  11,  1996 

Abstract 

The  transformation  of  potential  vorticity  within  and  stability  of  nonlinear  deep 
western  boundary  currents  in  an  idealized  tropical  ocean  are  studied  using  a  shallow- 
water  model. 

Observational  evidence  indicates  that  the  potential  vorticity  of  fluid  parcels  in 
deep  western  boundary  currents  must  change  sign  as  they  cross  the  equator,  but  this 
evidence  is  otherwise  unable  to  clarify  the  process.  A  series  of  numerical  experiments 
investigate  this  transformation  in  a  rectangular  basin  straddling  the  equator.  A  mass 
source  located  in  the  northwestern  corner  feeds  fluid  into  the  domain  where  it  is 
constrained  to  cross  the  equator  to  reach  a  distributed  mass  sink.  Dissipation  is 
included  as  momentum  diffusion.  The  Reynolds  number,  defined  as  the  ratio  of  the 
mass  source  per  unit  depth  to  the  viscosity,  determines  the  nature  of  the  flow,  and  a 
critical  value,  Rsc,  divides  its  possible  behavior  into  two  regimes.  For  Re  <  Rcc,  the 
flow  is  laminar  and  well  described  by  linear  theory.  For  Re  just  above  the  critical  value, 
the  flow  is  time-dependent,  with  cyclonic  eddies  forming  in  the  western  boundary 
current  near  the  equator.  For  still  larger  Reynolds  number,  eddies  of  both  signs 
emerge  and  form  a  complicated,  interacting  network  that  extends  into  the  basin  several 
deformation  radii  from  the  western  boundary,  as  well  as  north  and  south  of  the  equator. 

The  eddy  field  is  established  as  the  mechanism  for  potential  vorticity  transforma¬ 
tion  in  nonlinear  cross-equatorial  flow.  The  analysis  of  vorticity  fluxes  follows  from 
the  flux-conservative  form  of  the  absolute  vorticity  equation.  It  is  shown  that  the  zon- 
ally  integrated  meridional  flux  of  vorticity  across  the  equator  using  no  slip  boundary 
conditions  is  virtually  zero  even  in  the  strongly  nonlinear  limit  suggesting  that  the  ed¬ 
dies  are  extremely  efficient  vorticity  transfer  agents.  A  decomposition  of  the  vorticity 
fluxes  into  components  due  to  mean  advection,  eddy  transport,  and  friction,  reveals 
the  growth  with  Reynolds  number  of  a  turbulent  boundary  layer  that  exchanges  vor¬ 
ticity  between  the  inertial  portion  of  the  boundary  current  and  a  frictional  sub-layer 
where  modification  is  straightforward. 

A  linear  stability  analysis  of  the  shallow-water  system  in  the  tropical  ocean  exam- 
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ines  the  initial  formation  of  the  eddy  field.  The  formulation  assumes  that  the  basic 
state  is  purely  meridional  and  on  a  local  /-plane.  Realistic  western  boundary  current 
profiles  undergo  a  horizontal  shear  instability  that  is  partially  stabilized  by  viscosity. 
Calculations  at  several  latitudes  indicate  that  the  instability  is  enhanced  in  the  tropics 
where  the  internal  deformation  radius  is  a  maximum.  The  linear  stability  analysis  pre¬ 
dicts  a  length  scale  of  the  disturbance,  a  location  for  its  origin,  and  a  critical  Reynolds 
number  that  agree  well  with  numerical  results. 

Thesis  Supervisor:  John  Marshall  Title:  Professor,  Department  of  Earth,  Atmo¬ 
spheric,  and  Planetary  Sciences,  Massachusetts  Institute  of  Technology 
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Chapter  1 


Introduction 


During  intense  winter  storms  over  the  Labrador  and  Greenland  Seas,  water  cools  and 
convectively  adjusts,  sinking  to  great  depth.  In  so  doing,  it  begins  a  complicated 
journey  through  the  world’s  oceans,  where  it  eventually  rises  and  returns  to  these 
high  latitudes  only  to  begin  the  cycle  again.  Along  its  path,  most  of  this  water 
necessarily  crosses  the  equator,  at  least  twice  in  the  Atlantic,  and  perhaps  also  in 
the  Pacific  and  Indian  oceans.  Although  the  fact  of  deep  equatorial  crossing  is  well 
accepted,  the  details  of  the  process  remain  poorly  understood  from  both  observational 
and  dynamical  perspectives.  The  breakdown  of  the  geostrophic  approximation  at  low 
latitudes  leads  to  a  large  uncertainty  in  the  velocities  obtained  by  hydrographic  surveys 
in  the  vicinity,  and  the  strong  time-dependence  of  observed  currents  makes  mean  field 
calculations  from  current  measurements  particularly  questionable.  Analytical  studies 
have  successfully  modeled  near-equatorial  currents  under  very  restrictive  conditions, 
including  linear  and  purely  inviscid  dynamics.  However,  the  dynamical  regime  that 
most  likely  applies  to  the  ocean  is  the  intermediate  one,  in  which  inertial  motion  is 
significant  and  dissipation  weak  but  necessarily  non-zero. 

The  purpose  of  this  thesis  is  to  analyze  cross-equatorial  flow  under  these  more 
general  dynamical  conditions.  With  the  basic  assumption  that  fluid  must  cross  the 
equator,  the  question  of  interest  is  how.  The  problem  is  motivated  specifically  in  the 
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context  of  the  thermohaline  circulation.  The  surface  and  bottom  circulations  near  the 
equator  are  fundamentally  different  in  that  the  local  forcing  by  wind  and  topography 
contributes  to  the  vorticity  balance  and  thereby  changes  the  dynamical  constraints  of 
the  flow.  This  work  concentrates  on  non-locally  driven  currents,  with  the  deep  flow  of 
the  Atlantic  in  mind.  It  is  a  process-oriented  study,  numerical  in  nature  and  highly 
idealized  in  context,  intended  to  explain  dynamically  how  fluid  crosses  the  equator.  In 
particular,  it  determines  how  fluid  parcels,  carried  in  intense  cross-equatorial  currents 
and  constrained  by  only  weak  dissipation,  modify  their  potential  vorticity  along  their 
path.  The  thesis  will  show  that  this  process  takes  place  through  the  development  of  a 
turbulent  boundary  layer  that  resides  in  the  inertial  portion  of  the  western  boundary 
current  and  that  transfers  vorticity  from  this  outer  region  to  an  inner,  frictional  layer, 
where  it  is  removed  from  the  basin  through  lateral  dissipation. 

1.1  Motivation 

1.1.1  Hydrographic  Evidence  of  Cross-Equatorial  Flow 

A  considerable  body  of  evidence  supports  the  notion  that  fluid  travels  from  hemisphere 
to  hemisphere  in  the  ocean.  A  particularly  compelling  argument  is  found  in  any 
meridional  hydrographic  section  of  the  Atlantic  Ocean.  Shown  in  Figure  1.1  is  one 
such  section  (salinity)  from  the  Atlantic  GEOSECS  Atlas  {^Bainbridge  (1980)).  Using 
salinity  as  a  proxy  for  water  mass  reveals  that  several  different  water  masses  are 
found  in  the  Atlantic  both  north  and  south  of  the  equator.  In  particular,  the  high 
salinity  tongue  of  North  Atlantic  Deep  Water  (NADW)  plunges  to  depth  from  its  high 
latitude  formation  region  near  the  surface  southward  and  is  visible  between  depths  of 
1500  and  3500  meters  at  the  equator.  Sandwiching  this  NADW  are  two  layers  of  low 
salinity  fluid  of  southern  hemisphere  origin:  Antarctic  Bottom  Water  and  Antarctic 
Intermediate  Water. 

That  water  of  similar  origin  is  found  on  either  side  of  the  equator  invites  the 
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question  of  its  route,  which  can  not  be  inferred  from  this  section  alone,  but  which 
can  be  anticipated  from  dynamical  considerations  or  illustrated  through  additional 
observations. 

1.1.2  Dynamical  Motivation 

Viscous  limit 

Certainly  the  most  well-known  and  arguably  the  most  widely-accepted  theory  of  deep 
flow  in  the  ocean  is  that  of  Stommel  and  Arons  (1960a).  In  this  framework  as  in 
the  closely  related  homogenous  ocean  circulation  theories,  the  circulation  is  divided 
between  interior  and  boundary  layer  features.  The  interior  flow  is  governed  by  a  linear 
vorticity  balance, 

=  (1-1) 
where  f  is  the  Coriolis  parameter,  /?  its  meridional  derivative,  v  is  the  meridional  ve¬ 
locity,  H  the  fluid-layer  depth,  and  the  vertical  velocity  at  the  top  of  the  layer.  This 
upwelling  is  the  specified  forcing  of  the  system,  and  it  represents  a  distributed  vertical 
velocity  whose  basin-wide  integral  balances  the  volume  transport  of  fluid  sinking  in 
high-latitude,  convective  regions.  The  theory  assumes  a  level  bottom  so  the  vertical 
velocity  there  is  zero.  Appended  to  this  interior  flow  is  a  boundary  current  that  redis¬ 
tributes  mass  within  the  basin  and  which  can  be  shown  to  be  dynamically  consistent 
only  at  the  western  boundary. 

The  physics  of  this  boundary  current  was  left  unspecified  in  the  original  theory, 
but  its  existence  and  associated  transport  are  robust  features  of  the  model.  Assuming 
that  the  boundary  layer  dynamics  are  distinct  from  the  interior  flow,  the  boundary 
layer  transport  is  determined  entirely  by  the  interior  flow  and  spatial  distribution  of 
sources  and  sinks  that  supply  and  remove  fluid  in  the  deep  ocean.  The  boundary 
current  can  be  in  principle  either  inertial  or  frictional.  It  will  be  shown  below  that  a 
successful  model  of  cross-equatorial  flow  must  include  some  friction.  At  one  extreme  is 


the  purely  linear  example  of  the  Munk  (1950)  parameterization.  Modelling  the  oceanic 
dissipation  as  a  large-scale  diffusive  process,  the  boundary  layer  velocity  is  determined 

by 

ffv  =  (1.2) 


Here  meridional  motion  (i.e.,  planetary  vorticity  advection)  is  balanced  by  vorticity 
diffusion,  and  the  boundary  layer  has  a  zonal  scale, 

(1.3) 

The  dynamics  of  this  model  are  asymptotically  consistent  with  the  full  equations 
of  motion  provided  that  the  inertial  boundary  layer  scale, 

=  (1-4) 

where  U  is  the  scale  of  the  zonal  velocity,  is  small  compared  to  the  frictional  scale,  Sm- 
In  the  linear  theory,  currents  can  be  made  arbitrarily  weak  to  satisfy  this  constraint. 
At  the  equator  in  this  model,  a  weak  zonal  flow  in  the  interior  remains,  and  the 
boundary  layer  transport  reflects  the  geographical  location  of  the  sources  and  sinks, 
taking  fluid  from  hemisphere  to  hemisphere  as  required  by  the  global  balance.  There  is 
no  difficulty  modifying  vorticity  in  this  frictional  model.  One  conceptualization  of  the 
global  deep  circulation  of  the  ocean  is  shown  in  Figure  1.2  reflecting  these  dynamics. 
Notice  that  the  fluid  travels  unhindered  across  the  equator  in  each  ocean. 

Although  this  theory  is  successful  in  predicting  a  global  circulation,  it  requires  very 
weak  currents  to  validate  the  linear  dynamics  used.  From  the  full  vorticity  equation 
for  a  homogeneous  ocean, 

Ct+u-VC  +  /3v  +  if  +  C)^-u  =  AhV\  (1.5) 

where  (  is  the  relative  vorticity  of  the  flow,  it  is  clear  that  a  model  that  neglects 
relative  vorticity  requires  a  small  value  for  the  ratio, 

R  =  ^  (1.6a) 


(1.6b) 
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Figure  1.2:  Schematic  representation  of  the  global  abyssal  circulation  predicted  by  Stommel 

and  Arons.  Note  the  cross-equatorial  flow  indicated  in  each  ocean.  Adapted  from  Stommel  and 
Arons  (1960b). 


Here,  8^.  and  8y  represent  the  zonal  and  meridional  length  scales  over  which  the  fields 
vary.  It  is  possible  to  estimate  from  observations  these  scales  to  determine  if  the 
dynamical  assumptions  of  the  theory  are  consistent. 

Figure  1.3  presents  a  composite  velocity  profile  for  the  along-boundary  component 
in  the  tropical  Atlantic  ocean  at  1800  meters  depth  taken  from  SOFAR  floats  (Richard¬ 
son  and  Schmitz  (1993)).  There  exists  a  strong  southward  pecik  roughly  30  km  from 
the  continental  slope,  with  a  maximum  speed  of  26  cm/s.  In  their  text  the  authors  note 
that  the  full  data  were  averaged  into  10-km-wide  bins,  and  that  the  peak;  velocities 
reached  55  cm/s.  The  variation  of  this  velocity  scale  occurs  on  a  50  km  length  scale. 
The  velocity  perpendicular  to  the  boundary,  U,  and  the  meridional  length  scale,  Sy, 
are  not  well  known  from  observations,  but  using  V  =  .5  m/s,  8^^  =  50km,  ^  =  2x  10“^^^ 
111“^  s“^,  and  8y  =  1000  km,  gives 

R=l-  (1.7) 

This  0(1)  value  suggests  that  the  assumption  of  linear  dynamics  in  the  boundary 
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KILOMETERS  FROM  1800m 


Figure  1.3:  The  average  along-boundary  velocity,  transport  per  unit  depth,  and  number  of 

observations  from  the  1800  m  isobath  compiled  from  nine  near-equatorial  SOFAR  floats  at  1800 
meters  depth  in  the  western  Atlantic  Notice  the  intense  western  boundary  current  peaked  near 
a;  =  40  km  with  the  substantial  decline  in  velocity  immediately  adjacent  to  the  coast.  The  maximum 
velocity  in  the  DWBC  was  55  cm/s,  though  the  maximum  in  the  10  km-wide  averaged  bins  is  26 
cm/s.  Figure  adapted  from  Richardson  and  Schmitz  (1993). 
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current  is  inconsistent  with  available  observations.  Regardless  of  the  magnitude  of 
the  friction,  nonlinearity  must  be  considered  for  the  boundary  current  structure  and 
associated  dynamics  of  cross-equatorial  flow. 

Inviscid  limit 

Without  any  fundamental  understanding  of  the  magnitude  of  dissipation  in  the  ocean, 
it  is  tempting  to  take  an  opposing  view  and  consider  an  ocean  in  which  dissipation 
is  entirely  neglected.  Under  this  assumption,  the  central  dynamical  issue  of  cross- 
equatorial  flow  becomes  the  transformation  of  potential  vorticity.  This  quantity  is 
defined  by 

PV  = - (2f2  -f-  a;)  •  Vp,  (1-8) 

Po 

where  p  is  the  fluid  density,  po  a  reference  value,  f2  is  the  earth’s  rotation  vector,  and 
u  is  the  full  three-dimensional  curl  of  the  velocity  field.  The  potential  vorticity  is 
exactly  conserved  by  fluid  parcels  in  the  absence  of  forcing  or  dissipative  processes. 
As  a  result,  inviscid  and  unforced  meridional  motion  on  the  sphere  must  coincide  with 
changes  in  the  relative  vorticity  of  the  flow  or  in  the  stratification  of  the  fluid,  or 
both.  At  mid-latitudes  the  change  in  the  Coriolis  parameter  (/  =  2f2sin^,  where  <f>  is 
the  latitude)  is  relatively  small  compared  to  its  magnitude,  and  small  changes  in  the 
stratification  can  accomodate  large  meridional  excursions. 

However,  across  the  equator  the  projection  of  Q  on  the  local  vertical  changes 
sign,  and  the  fluid  must  compensate  by  creating  substantial  relative  vorticity  or  by 
inverting  the  density  field.  Far  from  the  western  boundary  current  and  far  from  the 
equator  such  as  in  the  large  portions  of  the  section  in  Figure  1.1,  relative  vorticity  is 
weak  compared  to  planetary  vorticity,  and  the  ocean  is  stably  stratified.  Assuming 
that  fluid  parcels  of  similar  origin  once  had  similar  potential  vorticity,  it  is  clear  that 
potential  vorticity  modification  is  fundamental  to  cross-equatorial  flow.  Somewhere 
along  their  route  from  hemisphere  to  hemisphere,  fluid  parcels  must  change  the  sign 
of  their  potential  vorticity.  As  this  change  can  only  be  achieved  through  some  form 
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of  forcing  or  dissipation,  the  process  of  cross-equatorial  flow  is  found  incompatible 
with  the  assumption  of  an  inviscid  ocean.  Ruling  out  local  forcing  since  this  study 
is  motivated  by  deep  ocean  dynamics,  dissipation  must  be  included  in  any  successful 
discussion  of  cross-equatorial  flow. 

As  shown  in  the  previous  section,  the  boundary  layer  dynamics  are  in  nature  at 
least  partially  inertial.  However,  some  dissipation  is  required  to  change  the  sign  of  the 
potential  vorticity.  Thus  the  ocean  operates  in  an  intermediate  regime,  both  inertial 
and  frictional. 


1.2  Background 

1.2.1  Observations 

Owing  to  the  impact  of  tropical  ocean  dynamics  on  global  weather  patterns,  there 
exists  a  sizeable  number  of  observations  of  the  tropical  ocean.  However,  these  stud¬ 
ies  are  primarily  devoted  to  the  surface  circulation  and  atmospheric  fluxes,  not  the 
deep  ocean  in  the  critical  western  boundary  region.  Some  observational  studies  have 
been  directed  toward  this  region  and  are  found  in  hydrographic  sections  and  float 
measurements.  Relying  on  thermal  wind  balance  to  deduce  currents  from  density  gra¬ 
dients,  indirect  estimates  from  hydrographic  surveys  are  limited  near  the  equator  by 
the  failure  of  the  geostrophic  approximation  there. 

Some  hydrographic  measurements  are  not  dependent  on  this  balance  and  there¬ 
fore  are  particularly  intriguing.  For  example,  the  measurement  of  CCI3F  on  the 
o'l.B  =  34.63kg  m“^  surface  in  the  the  core  of  Upper  North  Atlantic  Deep  Water 
(approximately  1600  meters)  by  Weiss  et  al.  (1985)  presents  an  indirect  measurement 
of  the  deep  equatorial  transport.  Reproduced  in  Figure  1.4  along  with  subsequent 
data  collected  five  years  later  (^Weiss  et  al.  (1989),  Weiss  et  al.  (1991),  and  Weiss  et 
al.  (1993)),  contours  of  high  chloroflouromethane  snake  down  the  western  boundary 
of  the  North  and  South  Atlantic  and  are  drawn  out  in  an  elongated  feature  directly  on 
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the  equator.  The  time-evolution  of  this  equatorial  tongue  indicated  by  a  comparison  of 
the  figures  is  unmistakable.  Clearly  some  fluid  in  the  deep  western  boundary  current 
is  directed  eastward  upon  reaching  the  equator,  while  the  remaining  fraction  travels 
directly  southward  across  the  equator. 

This  tracer  tongue  may  reflect  an  Eulerian  current  in  which  the  actual  circulation 
of  the  tropical  ocean  involves  an  along  equatorial  branch,  or  instead  it  may  indicate  a 
Lagrangian  pathway  from  the  boundary  to  the  interior  ocean  in  which  equatorial  waves 
draw  fluid  away  from  the  western  boundary  preferentially  at  the  equator.  The  original 
work  calculated  that  if  the  feature  reflected  the  former  possibility,  it  indicated  a  mean 
flow  for  the  boundary  current  of  approximately  1.4  cm  ( Weiss  et  ah  (1985)), 

McCartney  (1993)  compiled  data  from  3  hydrographic  surveys  that  bracket  the 
equatorial  Atlantic  in  an  effort  to  estimate  the  mean  circulation  of  the  deep  trop¬ 
ics.  Results  indicated  an  intense  deep  western  boundary  current,  transporting  35  x 
10  m  s  ^(355^1;)  of  NADW,  with  a  northward-flowing  counter-current  immediately 
east  of  the  DWBC,  returning  up  to  20  Sv,  From  property  distributions  and  geostrophic 
shear  calculations  at  6°N,  IFS,  and  between  2°  and  5°S  of  the  equator,  the  DWBC 
appeared  to  cross  the  equator  directly.  Although  the  counter-current  was  clearly  ob¬ 
served  in  both  hemispheres,  the  degree  of  its  connection  between  hemispheres  could 
not  be  established  from  the  data. 

SOFAR  floats  give  a  more  Lagrangian  impression  of  the  velocity  fleld  by  provid¬ 
ing  point  measurements  following  parcels  of  fluid  along  their  paths  {Richardson  and 
Schmitz  (1993)),  Shown  in  Figure  1.5  is  a  summary  of  displacement  vectors  from  float 
trajectories  observed  at  1800  meters  in  the  Atlantic.  It  is  this  data  set  that  also  corre¬ 
sponds  to  velocity  composition  in  Figure  1.3.  Multiple  trajectories  are  visible,  but  two 
characteristic  paths  stand  out.  One  set  of  trajectories  depicts  floats  travelling  directly 
across  the  equator,  at  least  to  about  5°S.  The  second  shows  an  along-equatorial  com¬ 
ponent,  with  floats  being  carried  both  east  and  west  up  to  1000  km  along  the  equator. 
Similar  to  the  tracer  data,  the  floats  suggest  a  bifurcation  of  fluid  at  the  equator. 
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CFC's  on  01.5  =  34.63 


Figure  1.4:  Composite  distribution  of  CCI3F  in  pMol  Kg“^  on  the  (T1.5  =  34.63  kg  m”^  surface 
(roughly  1600  meters  depth)  in  the  Atlantic  Ocean  determined  from  a)  TTO  cruises  (1981  and  1983) 
and  b)  five  years  later  from  SAVE  expedition  cruises  and  TTO  meaisurements  projected  forward  in 
time  based  on  their  apparent  CFG  ages.  Figure  courtesy  of  Ray  Weiss. 


Figure  1.5:  Displacement  vectors  of  SOFAR  float  trajectories  at  1800  meters  depth  between  Jan¬ 
uary  1989  and  October  1990.  Arrows  indicate  both  along-boundary  and  along-equatorward  motion. 
Figure  adapted  from  Richardson  and  Schmitz  (1993). 

The  variety  of  float  trajectories  in  Figure  1.5  also  indicates  considerable  variability 
in  the  region.  The  time-variability  of  the  tropics  has  been  noted  in  studies  of  the 
surface  circulation  particularly  in  the  Pacific,  but  such  behavior  is  characteristic  of 
the  tropical  ocean  at  all  depths.  Johns  et  al.  (1993)  found  60-70  day  fluctuations  of 
the  lower  NADW  boundary  current  properties  from  a  current  meter  moored  at  8°N, 
52°W.  Furthermore,  additional  SOFAR  floats  at  800  meters  demonstrated  loop-like 
behavior  indicating  sustained  eddy-like  features  propagating  in  a  northward-flowing 
intermediate-depth  western  boundary  current.  These  eddies  had  a  period  varying 
between  7  and  40  days. 

The  data  reflect  the  complicated,  time-dependent  circulation  of  the  region.  Even 
the  mean  path  of  fluid  parcels  in  the  tropical  circulation  is  not  well  understood  ob- 
servationally.  Although  it  seems  likely  that  the  fluid  travels  inter-hemispherically  in 
a  deep  western  boundary  current,  it  is  not  clear  whether  the  role  of  the  equatorial 
filament  is  significant.  In  addition,  it  is  not  evident  from  these  observations  how  fluid 
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modifies  its  potential  vorticity  along  its  path. 


1.2.2  Dynamical 

Several  dynamical  studies  of  equatorial  currents  have  contributed  to  the  general  back¬ 
ground  of  nonlinear  cross-equatorial  flow.  This  research  typically  falls  into  two  branches 
As  mentioned  above,  the  Stommel  and  Arons  (1960a)  theory  provided  the  initial  the¬ 
oretical  framework  for  considering  deep  flow  in  the  ocean,  although  it  applies  to  cur¬ 
rents  somewhat  weak  to  represent  oceanic  flows.  One  set  of  relevant  dynamical  studies 
emerges  from  this  formulation. 

Kawase  (1987)  extended  the  Stommel  and  Arons  work  by  parameterizing  the  large- 
scale  upwelling,  which  had  been  prescribed  in  the  earlier  theory,  so  that  the  vertical 
velocity  in  the  system  increased  with  layer  height  deviation  from  an  initial,  resting 
value.  The  upwelling  was  still  dependent  on  a  mixing  parameter.  A,  which  controlled 
the  sensitivity  of  the  upwelling  to  the  height  deviation.  In  so  doing,  this  model  ob¬ 
tained  a  steady-state  equatorial  tracer  tongue  for  large  values  of  A,  but  the  Stommel 
and  Arons  circulation  resulted  from  small  values  of  this  parameter.  Available  data 
suggested  that  the  value  required  for  a  tongue  corresponds  to  unrealistically  high  val¬ 
ues  of  the  vertical  mixing  coefficient  to  represent  features  of  low  vertical  mode.  In 
addition,  that  work  considered  the  time-dependent  behavior  of  the  spin-up.  It  showed 
that  an  along-equatorial  mass  flux  occurs  in  the  wake  of  an  equatorial  Kelvin  wave 
resulting  from  a  change  in  the  source  strength.  Kawase  concluded  that  the  equatorial 
signal  must  result  from  time-dependent  features  of  the  circulation  and  not  from  a 
steady-state  balance.  The  dynamics  of  this  model  were  primarily  linear. 

A  subsequent  contribution  by  Kawase  et  al.  (1992)  showed  the  steady  behavior 
of  a  full  three-dimensional  general  circulation  model  for  this  problem.  Forced  by  an 
imposed  temperature  anomaly  at  the  northern  boundary,  the  authors  found  again 
that  the  final  state  of  the  system  consisted  predominantly  of  a  time-averaged,  cross- 
equatorial  flow.  Able  to  extend  the  system  into  a  more  nonlinear  range  than  formerly. 
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the  study  also  demonstrated  the  development  of  time-dependent,  eddy-like  features  in 
the  boundary  region. 

Springer  and  Kawase  (1993)  explored  the  connection  between  a  nonlinear  western 
boundary  and  equatorial  currents  in  an  experiment  designed  to  prevent  the  full  spin- 
up  of  the  deep  equatorial  ocean.  A  sponge  layer  was  placed  at  the  eastern  boundary  to 
accept  the  initial  Kelvin  wave  involved  in  the  spin-up  of  the  tropical  ocean  and  prevent 
the  mid-latitude  interior  from  reaching  the  Stommel  and  Arons  steady  state.  Thus 
the  model  reached  a  quasi-steady  state  that  included  a  western  boundary  current  in 
one  hemisphere  coupling  to  an  along-equatorial  meandering  flow. 

The  second  branch  of  relevant  dynamical  studies  emerges  from  models  of  the  sub¬ 
tropical  gyre,  and  originates  with  the  numerical  demonstration  by  Bryan  (1963)  that 
the  weakly  nonlinear  Munk  boundary  layer  is  unstable.  More  recent  studies  (e.g., 
lerley  (1987),  Kamenkovich  et  al.  (1995),  and  Sheremet  et  al.  (  1995))  have  extended 
the  study  to  include  very  nonlinear  flows.  These  authors  also  find  a  break-down  of 
the  steady  laminar  Munk  boundary  layer  and  for  particularly  strong  nonlinearity,  the 
development  of  an  inertial  recirculation  in  one  region  of  the  basin.  This  recirculation 
is  a  persistent  feature  in  the  mid-latitude  gyre,  necessary  to  enhance  the  transport  in 
the  boundary  layer  to  dissipate  vorticity  in  the  region.  Although  it  is  not  necessary 
to  change  the  sign  of  the  potential  vorticity  in  a  mid-latitude  model,  it  is  necessary 
to  balance  vorticity  input  by  the  imposed  wind-stress  and  reach  a  statistically  steady 
state. 

The  instability  of  the  Munk  boundary  layer  at  mid-latitudes  has  been  examined 
through  a  linear  stability  analysis  by  lerley  and  Young  (1991)  and  Cessi  and  ler¬ 
ley  (1993).  These  studies  demonstrate  that  the  boundary  layer  is  unstable  to  small 
perturbations  at  Reynolds  numbers  between  20  and  50,  depending  on  the  boundary 
conditions  and  model  configuration.  They  define  the  instability  as  viscous  because  the 
trapped,  growing  modes  result  from  the  viscous  term  in  their  linear  stability  equation. 
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1.3  Outline  of  the  Thesis 


This  thesis  examines  the  dynamics  of  western  boundary  currents  constrained  to  cross 
the  equator  and  is  divided  in  principle  into  two  portions.  The  first  portion  focusses  on 
the  central  question  of  the  thesis:  how  is  the  potential  vorticity  modified  in  a  strongly 
nonlinear  western  boundary  current  to  enable  cross-equatorial  flow?  Chapters  2,  3, 
and  4  directly  address  this  problem  through  the  development  of  a  numerical  model 
and  its  diagnosis  that  shows  that  this  process  is  carried  out  through  the  exchange  of 
vorticity  between  frictional  and  turbulent  boundary  layers.  The  focus  of  Chapter  5 
is  more  indirectly  related  to  this  question,  discussing  the  cause  of  this  dynamically 
essential  eddy  field. 

Chapter  2  presents  an  overview  of  the  model  used  in  this  thesis.  It  begins  with 
a  scaling  of  contributions  to  Ertel’s  potential  vorticity  which  determines  that  the 
shallow-water  model  is  the  simplest  representation  of  the  full  Navier-Stokes  equations 
that  adequately  models  cross-equatorial  flow.  A  subsequent  description  outlines  the 
shallow-water  model  itself  and  the  configuration  of  the  domain  and  forcing  used  in 
this  study.  The  model  is  highly  idealized  with  a  rectangular  domain  and  flat  bot¬ 
tom,  and  it  succesfully  constrains  fluid  to  cross  the  equator  through  a  mass  balance 
between  a  source  in  the  northern  hemisphere  and  a  sink  that  is  either  distributed 
uniformly  around  the  basin,  or  localized  in  the  southern  hemisphere.  One  critical 
non-dimensional  parameter  that  determines  the  overall  character  of  the  system  is  the 
Reynolds  number,  which  is  defined  a^  the  ratio  of  the  overall  mass-flux  through  the 
system  per  unit  depth  to  the  viscosity.  Finally,  the  numerical  model  is  introduced 
with  a  brief  description  of  its  properties  and  necessary  boundary  conditions. 

The  numerical  experiments  for  a  range  of  Reynolds  numbers  are  discussed  in  Chap¬ 
ter  3.  At  small  Re,  the  system  obeys  linear  dynamics  and  the  connection  to  the 
Stommel  and  Arons  (1960a)  work  is  established.  At  a  critical  Re,  the  system  be¬ 
comes  unstable  with  eddies  developing  immediately  north  of  the  equator,  propagating 
southward  and  ultimately  dissipating  in  the  southern  hemisphere.  The  sense  of  the 
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circulation  of  these  eddies  follows  from  the  vorticity  of  the  fluid  at  the  mass  source. 
For  still  higher  Reynolds  numbers,  the  basin  is  filled  with  time-dependent  features  of 
both  signs.  Although  instantaneous  pictures  of  the  height  and  velocity  fields  show 
a  complicated  pattern  of  interacting  vortices,  the  time-mean  picture  reveals  beneath 
the  eddy  field  a  well-defined  western  boundary  current  transporting  fluid  across  the 
equator. 

The  diagnosis  of  these  experiments  is  the  subject  of  Chapter  4.  After  introduc¬ 
ing  the  J-vector  as  the  primary  tool  to  investigate  the  dynamics,  examples  of  its 
application  to  the  simple  linear  and  weakly  nonlinear  systems  illustrate  analytically 
how  vorticity  fluxes  through  a  basin  strongly  constrained  by  friction.  Analysis  of  the 
strongly  nonlinear  experiments  follows  and  shows  the  growth  with  Reynolds  number 
of  a  turbulent  boundary  layer  that  transports  vorticity  between  the  inertial  portions  of 
the  full  boundary  current  and  a  frictional  sublayer  adjacent  to  the  boundary  where  its 
transformation  is  clear.  Thus  the  eddies  have  the  essential  dynamical  role  as  the  agents 
of  vorticity  modification  that  allows  cross-equatorial  flow.  The  chapter  also  dissects 
the  time-mean  boundary  layer  structure  and  includes  a  discussion  of  the  sensitivity  of 
the  results  to  different  boundary  conditions  and  model  configurations. 

Chapter  5  examines  the  cause  of  the  eddy  field  that  emerges  in  the  numerical  model 
runs  through  a  linear  stability  analysis  of  the  western  boundary  current.  Specifically 
three  possible  mechanisms  are  investigated,  inviscid  shear  instability,  viscous  instabil¬ 
ity,  and  inertial  instability.  The  formulation  of  the  instability  problem  is  followed  by 
several  sections  on  the  instability  under  a  range  of  conditions.  Section  5.3.3,  which 
considers  geostrophic  perturbations,  examines  the  most  idealized  model  and  provides 
the  simplest  mathematical  and  physical  description  of  the  instability,  finding  that  the 
western  boundary  current  suffers  a  shear  instability  resulting  from  the  horizontal  varia¬ 
tion  in  the  potential  vorticty  structure.  This  section  also  illustrates  that  the  instability 
is  tropically  enhanced  because  the  maximum  in  the  deformation  radius  occurs  at  the 
equator. 
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Chapter  2 


The  Model 

2.1  Scaling  Ertel’s  Potential  Vorticity 

The  dynamical  problem  of  inertio-frictional  cross-equatorial  flow  centers  upon  the  bal¬ 
ance  between  conservation  and  modification  of  potential  vorticity.  The  dynamics  of 
this  flow  could  be  studied  with  a  range  of  models  of  varying  complexity.  However, 
as  in  other  types  of  analysis,  it  is  desirable  to  condense  the  necessary  features  of  the 
full  system  into  a  simple  form  that  focuses  attention  on  those  fundamental  features. 
To  determine  the  simplest  model  that  contains  the  necessary  physics  to  study  cross- 
equatorial  flow,  a  comparison  of  the  magnitude  of  different  terms  of  the  total  potential 
vorticity  highlights  the  components  of  the  full  momentum  equations  that  govern  in- 
viscid  flow.  ErtePs  potential  vorticity,  Equation  1.8,  can  be  re-expressed  in  Cartesian 
coordinates  and  using  subscripts  to  denote  partial  derivatives, 

PV  —  --  {{Wy  -  •U^)/9x  +  (Wz  -Wx+  2fi)  Py  +  (v^  -Uy+  fiv)  P z}  -  (2.l) 

P 

The  equatorial  /0-plane  approximation  in  which  f3  =  is  the  local  gradient  of 

planetary  vorticity  at  the  equator,  t/  =  0,  has  been  used.  The  zonal,  meridional,  and 
vertical  components  of  velocity  are  represented  by  u,  u,  and  u;,  respectively.  Although 
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limited  observations  relating  to  cross-equatorial  flow  exist,  it  is  possible  to  combine 
different  measurements  to  estimate  reasonable  scales  for  the  various  components  of 
Equation  2.1.  As  in  the  introduction,  I  approximate  the  along-boundary  component 
of  the  flow  to  be  meridional  and  take  from  Figure  1.3  F  =  50cm  s~^  and  =  100  km. 
I  assume  that  the  vertical  scale  is  D  =  500  meters  and  that  the  vertical  velocity  scales 
according  to  continuity:  W  <  Since  the  mean  zonal  velocity  and  meridional 

scale  are  not  well  determined  by  observations,  I  assume  an  upper  limit  on  their  values: 
U  =  V  and  6y  —  8x-  Density  variations  can  be  estimated  to  be  py  =  5  x  10~^A:p/m^ 
and  px  —  3  X  10^^ kg fm'^  {McCartney  (1993)).  Given  these  values,  all  components 
of  relative  vorticity  are  approximately  of  the  same  scale,  but  the  contributions  of 
the  horizontal  components  to  the  potential  vorticity  are  small  compared  to  that  in 
the  vertical.  Furthermore,  using  2f2(=  1.4  x  10“^s“^)  and  /3  =  2  x 
the  contribution  by  the  horizontal  component  of  planetary  vorticity  is  also  small  in 
Equation  2.1  whereas  that  due  to  the  vertical  component  can  not  be  neglected. 

Equation  2.1  then  reduces  to 

PV  =  — -  {{(3y  +  Vx  -  Uy)  px}  (2.2) 

Po 

which  corresponds  directly  to  the  shallow-water  model.  A  single  fluid  layer  contains 
the  necessary  physics  to  model  observed  cross-equatorial  flows,  neglecting  interactions 
between  layers.  It  is  the  difference  between  the  internal  compensation  for  changes 
in  planetary  vorticity  through  stretching  versus  that  through  changes  in  relative  vor¬ 
ticity  that  makes  the  equator  unique  among  latitudes,  and  the  shallow-water  model 
comprises  both  of  these  elements. 

It  is  reasonable  to  inquire  under  what  conditions  the  above  scaling  breaks  down. 
For  flows  with  strong  meridional  density  variations,  in  a  narrow  band  about  the  equa¬ 
tor,  the  horizontal  component  of  vorticity  must  be  included  in  the  dynamics.  Naturally 
the  shallow-water  model  is  also  insufficient  to  resolve  flows  with  large  vertical  veloci¬ 
ties.  The  use  of  the  present  model  omits  the  possibility  of  density  inversion. 
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2.2  The  Shallow- Water  Model 


The  shallow-water  model  is  a  simplification  of  the  full  Navier-Stokes  equations  that 
describes  the  flow  of  a  single  layer  of  fluid,  the  active  layer,  generally  bounded  by 
a  motionless  layer  of  differing  density,  either  above  or  below.  Opposite  the  inactive 
layer  is  usually  a  solid  boundary  that  contains  the  fluid  with  walls  bounding  the  basin 
laterally.  For  a  system  in  which  vertical  scales  are  small  compared  to  horizontal  scales 
and  density  variations  within  the  fluid  can  be  approximated  to  a  narrow  interface,  the 
primitive  equations  of  motion  and  continuity  simplify  to 


Ut -f  u  •  Vu -|- /k  X  u  =  —giVh  +  F  (2.3a) 

ht  +  u-Vh  +  hV  -u  =  -iw*.  (2.3b) 

Here,  u  =  (u,v)  is  the  horizontal  velocity  vector,  gt  =  is  the  gravitational  acceler¬ 
ation,  reduced  by  the  ratio  of  density  variations  at  the  interface  to  the  density  of  the 
layer  itself,  and  h  is  the  depth  of  the  fluid.  The  forcing  and  dissipation  are  included  in 
body  force  terms,  F,  and  in  tn*  which  represents  the  cross-isopycnal  mass  flux  between 
the  active  and  inactive  layers.  The  system  as  described  in  this  thesis  is  laid  out  in 
Cartesian  coordinates. 

If  the  flow  is  strongly  constrained  by  the  Coriolis  acceleration,  velocities  are  ap¬ 
proximately  depth-independent.  One  measure  of  this  strength  is  the  Rossby  number, 
Ro  =  which  compares  advective  to  planetary  time  scales.  Mid-latitude  currents 
away  from  the  western  boundary  typically  are  characterized  by  a  small  Rossby  num¬ 
ber.  This  study,  however,  considers  strong  flows  at  low  latitudes,  and  therefore  the 
Rossby  number  is  0(1).  However,  strong  currents  are  still  in  cross-stream  geostrophic 
balance  even  near  the  equator.  Any  vertical  structure  in  the  velocity  field  is  neglected 
by  the  shallow- water  equations. 

The  above  equations  relate  to  a  homogeneous  layer.  The  model  is  general  in  the 
sense  that  it  can  represent  the  full  ocean  as  a  single  layer,  in  which  case  the  bounding. 
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inactive  fluid  is  air,  and  gt  reduces  to  g.  In  this  barotropic  limit,  layer  depth  variations 
are  typically  quite  small  since  even  small  variations  produce  large  pressure  gradients 
and  correspondingly  large  velocities.  For  cases  in  which  the  interface  is  internal  as  in 
this  study,  g'  <  g,  and  the  system  is  considered  a  l|  layer  model.  The  shallow- water 
model  is  also  general  in  that  it  can  represent  an  active  layer  sandwiched  between  two 
quiescent  fields,  as  in  Pedlosky  (1987b).  It  is  this  configuration  that  best  depicts  the 
deep  flows,  far  above  the  ocean  bottom. 

In  nature,  the  ocean  is  driven  by  atmospheric  forcing,  by  winds  and  by  thermal 
and  salinity  variations,  and  dissipation  is  required  to  balance  this  input.  The  eff^ect  of 
winds  is  included  in  a  shallow-water  model  through  a  wind-stress,  r,  in  the  momentum 
equations.  Thermal  and  saline  effects  can  be  approximated  by  a  single  buoyancy 
forcing  in  the  density  equation  of  a  full  three-dimensional,  baroclinic  model.  Through 
atmospheric  forcing  and  small-scale  mixing,  fluid  density  is  transformed  from  one 
value  to  another,  resulting  effectively  in  a  cross-isopycnal  mass  flux.  This  influence  is 
incorporated  into  the  shallow-water  system  through  a  term,  B,  in  the  height  equation. 
Most  generally,  dissipation  is  included  in  both  the  momentum  equations,  representing 
a  large-scale  viscosity,  and  in  the  height  equation,  as  due  to  diffusion  of  density  by 
vertical  mixing. 

With  these  additions,  the  equations  take  the  form 


Ut -f- u  •  Vu -|- /k  X  u  =  —gfVh  +  T  —  vise  (2.4a) 

h't  -k  u  •  Vh  -f-  /iV  •  u  =  B  —  Di  f  f .  (2.4b) 

Since  the  present  study  is  directed  toward  deep  flows  in  the  ocean,  the  wind-stress 
is  zero,  and  the  vertical  mass  flux,  B,  is  assumed  to  correspond  to  high-latitude  forcing 
in  regions  where  convection  occurs.  Thus,  let 

B  =  B{x,y)  (2.5) 
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where  this  function  is  assumed  to  be  zero  except  in  a  small  region  near  the  northern¬ 
most  latitude  of  the  model.  The  function  is  constrained  by  a  mass  balance  such  that 

JhdA  =  j  jBdA-j  j  DiffdA  (2.6) 

which  must  be  zero  in  the  steady  limit. 

The  dissipation  in  the  real  ocean  is  poorly  understood,  and  several  research  groups 
are  active  in  its  study,  from  both  observational  and  modeling  perspectives.  Rather 
than  focussing  on  the  intricacies  that  undoubtedly  develop  with  different  forms  of 
dissipation  and  determining  which  form  represents  the  most  realistic,  this  study  in¬ 
vestigates  only  the  simplest  of  parameterizations,  with  general  conclusions  drawn  that 
may  be  useful  when  the  subject  of  mixing  is  better  understood.  Viscosity  is  included 
as  lateral  friction  .  Specifically, 

vise  =  AhV^u,  (2.7) 

and  diffusion  is  a  specified  function  of  position, 

Diff  =  Diff(x,y).  (2.8) 

Observations  suggest  that  there  is  substantial  spatial  structure  in  the  vertical  mixing 
of  the  ocean.  This  structure  seems  to  depend  on  the  proximity  of  topographic  features 
which  enhance  the  internal  wave  field  in  its  vicinity  {J.  Toole  (personal  communication) 
and  Kunze  and  Sanford  (1996)).  More  than  likely,  it  also  depends  on  the  internal 
dynamics  of  the  system,  such  as  the  strength  of  the  flow  field  and  mass  distribution. 
However,  this  dependence  is  difficult  to  estimate  with  any  confidence.  As  a  result,  the 
mixing  is  assumed  to  be  either  uniform  or  localized.  With  buoyancy  and  mixing  pre¬ 
determined  functions  of  location,  the  net  mass  flux  in  the  model  is  specified  entirely. 
However,  the  pathways  by  which  fluid  moves  from  the  source  at  high  latitudes  to  the 
sink  is  internally  determined  and  this  is  the  principle  focus  of  the  thesis. 

No  normal  mass  flux  (no  flow)  conditions  are  applied  at  the  lateral  boundaries  in 
all  runs.  No  slip  is  used  in  most  runs,  with  some  attention  paid  briefly  to  free  slip. 
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Thus,  on  each  boundary, 

u  •  n  =  0 

(2.9) 

and,  either 

o 

II 

(2.10) 

or 

Vu  •  n  =  0 

(2.11) 

where  n  and  i  are  unit  vectors  normal  and  tangent  to  the  boundary. 


2.3  Model  Configuration  and  Scaling 

To  study  cross-equatorial  flow,  I  conflgure  a  maximally  simple  geometry,  a  rectangular 
basin  with  vertical  lateral  boundaries,  straddling  the  equator,  on  the  equatorial  /3- 
plane.  The  basin  extends  a  distance  zonally,  and  Ly  meridionally,  with  y  =  0 
corresponding  to  the  central  latitude  as  in  Figure  2.1.  Although  the  shallow- water 
equations  can  not  distinguish  between  mass-fluxes  through  either  an  upper  or  lower 
interface  ^  the  dynamics  of  the  vorticity  equation  respond  only  to  vortex  stretching 
not  its  location  —  the  inflow  has  meaning  with  regard  to  the  sign  convention  of  the 
net  mass  flux,  w^.  As  Equations  2.3a— 2.3b  are  written  and  assuming  the  active  layer 
rests  beneath  the  inactive  layer,  a  positive  corresponds  to  upwelling  through  the 
interface,  and  a  negative  value  implies  downwelling. 

In  this  configuration,  it  is  convenient  to  non-dimensionalize  so  that  velocities  are 
scaled  by  the  gravity  wave  speed,  U  V  =  Cq  =  ^/^,  horizontal  lengths  by  the 
deformation  radius,  =  Sy  =  Ld  =  ^,  W  =  Ah  =  and  the  time- 

scale,  T  =  Dimensional  units  are  denoted  with  hats.  With  these  deflnitions,  the 
shallow  water  equations  take  the  simple  form 

ut -I- u  •  Vu -fyk  X  u  =  -S/h  +  AnV^u 
-f-  u  •  V/i  -f-  /iV  •  u  =  — Wqu;*. 
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(2.12a) 

(2.12b) 


Lx 


Figure  2.1:  Schematic  of  the  model  under  conditions  of  strong  inflow  in  the  northwest  corner  (x) 
and  weak  uniform  upwelling  throughout  the  basin  (shaded). 

Wo  =  represents  the  basin-wide  average  strength  of  the  upwelling.  Stated 

differently,  So  is  the  sourqe  volume  that  is  forced  through  the  system,  injected  into 
the  northwest  corner,  and  removed  elsewhere  in  the  basin.  Lj:  and  Ly  are  the  non- 
dimensional  basin  dimensions.  The  volume  flux.  So,  is  scaled  in  units  of  HcLd,  For 
H  =  500m,  c  =  Im/s,  and  Lo  =  250A:m,  =  0.2  HcLj)  corresponds  to  an  actual 

transport  of  25  Sv. 

With  this  non-dimensionalization,  the  four  external  parameters  of  the  original 
system,  /?,  g\  Ah,  and  Wo,  conveniently  reduce  to  two,  the  forcing  and  dissipation. 
The  range  of  model  behavior  is  more  easily  spanned. 

Introducing  the  the  top-hat  function. 


1  if  la:]  <  I 
0  if  |x|  > 


(2.13) 


35 


the  forcing  used  in  most  experiments  of  the  thesis  hcis  the  form 


= 


Lx  Ly 


AxAy  V  Aa: 
Additional  experiments  use  more  localized  upwelling, 


(2.14) 


Parameters,  xO,  xl,  etc.  shift  the  central  location  of  the  source  and  sink,  and  the 
dimensions  Ax  and  Ay  give  the  widths  of  those  regions.  Referring  to  experiments 
with  the  former  mass-flux  distribution  as  the  “uniform”  upwelling  case  is  a  partial 
misnomer  since  the  mass  source  is  also  included  in  to*.  However,  since  the  interior  flow 
is  determined  by  the  local  vorticity  balance  (i.e.,  the  local  upwelling),  it  is  convenient 
hereafter  to  distinguish  the  two  forcing  distributions  as  “uniform”  (Equation  2.14)  and 
non-uniform”  or  “localized”  (Equation  2.15).  An  example  of  the  localized  upwelling 
configuration  is  given  in  Section  4.2,  where  this  distribution  is  applied. 


2,4  The  Vorticity  Equation  and  Reynolds  Num¬ 
ber 

Although  two  parameters,  Wq  and  Ah,  govern  Equations  2.12a  -  2.12b,  the  dynamics 
may  be  simplified  further  by  considering  a  single  parameter.  To  this  end  it  is  useful 
to  derive  the  vorticity  equation,  obtained  by  taking  the  curl  of  the  shallow-water 
momentum  equations: 

— -I-U- VC  +  v-f-(C  +  2/)V  u=:  (2.16) 

As  in  the  introduction,  boundary  layer  scales  can  be  estimated  assuming  different 
steady-state  balances  in  Equation  2.16.  For  example,  a  balance  of  planetary  vorticity 
advection  and  vorticity  diffusion  gives  the  Munk  boundary  layer  scale, 

Sm  =  Aff.  (2.17) 
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The  balance  of  relative  and  planetary  vorticity  advection  leads  to  the  inertial  boundary 
layer  scale, 

6i  =  VU  (2.18) 

where  U  is  the  scale  of  the  zonal  velocity  within  the  boundary  layer.  The  forcing 
parameter,  Sq,  relates  to  the  inertial  boundary  layer  scale  by  setting  the  magnitude 
of  the  zonal  current  to  which  the  layer  matches.  Under  conditions  of  uniform  forcing, 
for  example. 

As  a  result,  the  two  parameter  system  of  forcing  and  dissipation  can  be  restated  in 
terms  of  Munk  and  inertial  boundary  layer  scales. 

An  additional  length  scale  is  derived  assuming  a  balance  of  nonlinearity  and  friction 
in  Equation  2,16.  This  inertio-frictional  scale  is  defined 


^  —  Y  (2.20) 

This  balance  arises  when  the  full  boundary  layer  structure  is  predominantly  inertial. 
Even  for  a  very  small  viscosity,  friction  must  ultimately  become  significant  in  the 
balance  at  some  scale,  though  one  smaller  than  the  Munk  scale. 

It  is  convenient  to  define  the  Reynolds  number,  which  represents  the  relative 
strength  of  nonlinearity  and  dissipation.  The  ratio  of  these  two  terms  in  Equation  2.16 


gives 


Rci  —  -- — 

Ah 


(2.21) 


where  S  represents  the  zonal  scale  over  which  vorticity  variations  occur.  Taking  8  =  8i 
as  for  example  in  Kamenkovich  et  al.  (1995)  gives 

Re.=  0.  (2.22) 

Although  this  quantity  is  related  to  the  external  parameters  of  the  system  through 
Equations  2.17  and  2.19,  the  definition  of  the  inertial  boundary  layer  width  remains 
undesirably  dependent  on  the  distribution  of  the  upwelling  in  the  interior. 
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A  second  measure  of  relative  strengths  of  nonlinearity  and  dissipation  can  be  de¬ 
fined  as  the  ratio  of  the  source  strength  per  unit  depth  to  the  viscosity: 

50 


Reo  = 


AhH' 


(2.23) 


I  prefer  the  second  measure,  which  follows  lerley  and  Young  (1991),  since  it  is  unam¬ 
biguously  defined  by  the  input  parameters.  Future  references  to  the  Reynolds  number 
will  drop  the  subscript,  2. 


2.5  The  Numerical  Model 

2.5.1  Spatial  Discretization 

To  study  the  system  numerically.  Equations  2.12a  and  2.12b  must  be  discretized,  and 
several  choices  are  available.  I  employ  a  modified  version  of  the  potential  enstrophy 
conserving  model  of  Sadoumy  (1975)  for  this  study.  In  this  model  a  numerical  form 
of  the  mass-weighted  square  of  the  potential  vorticity  is  approximately  conserved.  A 
comparison  to  a  similar  energy  conserving  model  showed  that  in  the  inviscid  limit  a 
numerical  scheme  that  does  not  conserve  potential  enstrophy  results  in  the  unphysical 
growth  of  rotational  energy  at  small  scales,  leading  to  unrealistic  flows  and  unre¬ 
alistic  energy  dissipation.  Potential  enstrophy  conservation  avoids  this  small  scale 
catastrophe.  Another  useful  model  is  the  potential  enstrophy  and  energy  conserv¬ 
ing  finite-difference  scheme  of  Arakawa  and  Lamb  (1978),  described  in  Haltiner  and 
Williams  (1980).  With  more  conserved  quantities,  this  scheme  potentially  produces 
flows  that  are  more  consistent  with  the  continuous  system.  However,  the  model  is 
more  complicated  and  therefore  more  difficult  to  diagnose  accurately.  The  Sadourny 
model  has  the  advantage  that  it  is  relatively  simple  and  slightly  faster  to  run.  I  have 
examined  both  models  and  qualitative  comparisons  of  the  output  show  good  agree¬ 
ment  in  the  gross  behavior  of  the  system.  The  inclusion  of  dissipation  in  this  work, 
even  though  quite  weak,  renders  the  conservation  properties  of  the  model  less  vital 
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A  O  A  O  A  O  A 


Figure  2.2:  Schematic  of  the  grid  used  in  the  model, 
than  in  inviscid  studies. 

The  numerical  model  is  laid  out  on  a  C-grid,  on  which  dependent  variables  are 
staggered.  Consider  the  basin  divided  up  into  a  series  of  cells  of  uniform  dimension, 
as  in  Figure  2.2  with  layer  depth,  /i,  at  their  centers.  Vorticity  and  potential  vorticity 
lie  at  cell  corners,  and  velocities,  u  and  v,  are  directed  perpendicular  to  and  located 
at  the  midpoint  of  the  bounding  line  segments  of  each  cell. 

Notationally,  varibles  axe  defined  where  i  repre¬ 

sents  the  zonal  index  and  j,  the  meridional.  It  is  convenient  to  create  averaging  and 
differencing  operators,  defined  at  midpoints  of  connecting  line  segments,  and  a  Lapla- 
cian  operator,  defined  at  the  point  at  which  it  operates.  For  example,  in  the  case  of 
relative  vorticity,  let 

C  i+^.j  =  2 
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(2.24a) 

(2.24b) 


and  similarly  in  y,  and 


^  Ci,j  —  (^xCt,j)  d"  •  (2.24c) 

Ax  is  the  grid  spacing  in  the  x-direction.  It  can  be  shown  the  differential  and  averaging 
operators  commute  with  themselves,  so  that  multiple  operations  may  be  taken  in  any 
order. 

Using  this  notation,  the  discretization  takes  the  form 


(2.25a) 

=  ^hA'u 

(2.25b) 

ht  +  SxU  +  6yV  —  — iw*. 

(2.25c) 

where  the  potential  vorticity  is  defined 

f  +  C 

r» 

(2.26) 

and 

(  =  SxV  —  8yU 

(2.27a) 

II 

i 

(2.27b) 

V  =  Vv 

(2.27c) 

<t>  =  +  ^  . 

(2.27d) 

As  in  Adcroft  (1995),  several  rules  can  be  determined  from  the  operators  2.24a- 
2.24c.  Letting  k  represent  either  horizontal  direction,  and  variables  'tp  and  rj  be  general, 

(2.28a) 

4(V’  v)  =  i’^kV  +  rfSkip  (2.28b) 

and  the  square  of  the  geometric  mean  of  a  variable  is 

=  V’fc^fc+i- 


(2.28c) 


The  potential  enstrophy  of  the  continuous  system  is  defined,  ^  =  hq^,  and  the  following 
equation  can  be  derived  from  the  discretization  2.25a-2.25c  with  Ah  =  =  0, 

^  {q^TT)  +  S,  {iq  tl  r")  +  S,  (tg  T’)  =  0,  (2.29) 

Thus  an  approximate  form  of  the  potential  enstrophy  is  conserved. 

Details  of  the  numerical  model  time-stepping  and  grid-stretching  can  be  found  in 
Appendix  A. 


2.5.2  Numerical  Boundary  Condition 


The  boundary  conditions  for  the  continuous  system  consist  of  no  normal  flow  through 
the  solid  boundaries,  and  no  slip  or  free  slip  at  the  batsin  walls.  However,  the  discretiza¬ 
tion,  Equations  2.25a  -  2.25c,  requires  the  specification  of  the  potential  vorticity  at  all 
grid  points  on  the  boundary  as  well  as -within  the  domain,  and  therefore  an  additional 
boundary  condition  is  required  to  step  forward  the  model.  Although  the  no  flow, 
no/free  slip  conditions  determine  the  absolute  vorticity  on  the  boundary,  they  do  not 
determine  the  layer  height.  This  must  be  determined  separately,  but  it  is  desired  to 
be  consistent  with  the  continuous  system. 

Consider  the  zonal  momentum  equation, 

Du  ^  idh  ^ 

=  (2.30) 

where  ^  represents  the  material  derivative.  Setting  ti  =  0  as  is  true  at  meridional 
boundaries  and  rearranging  terms  leaves 


dh  1  /  ,  d^u\ 


(2.31) 


The  meridional  velocity  is  determined  by  the  no/free  slip  condition,  but  the  curva¬ 
ture  of  the  zonal  velocity  is  not.  I  choose  the  simplest,  though  somewhat  arbitrary, 
boundary  condition  that  sets  the  gradient  of  the  curvature  to  zero,  namely 

d^u 


dx^ 


=  0 


(2.32) 


41 


at  the  meridional  boundaries,  and  similarly  in  v  at  the  zonal  boundaries. 

A  second  boundary  condition  that  may  be  viable  would  be  to  step  forward  the 
height  equation  at  the  boundary, 

Dh 

—  +  W  •  u  =  -Wow^.  (2.33) 

Although  I  have  not  implemented  this  form  for  the  boundary  condition,  I  can  assess 
one  measure  of  consistency  between  the  chosen  condition  with  the  one  that  is  exactly 
consistent  with  the  continuous  system.  If  the  two  conditions  were  inconsistent,  then 
the  model  would  not  conserve  mass.  This  represents  an  integrated,  rather  than  local 
check,  but  I  have  compared  mass  starting  an  integration  and  mass  after  3000  time 
units,  and  found  the  two  to  match  exactly. 

2.5.3  Parallel  Programming 

Although  the  model  is  general  and  has  run  on  various  serial  platforms,  it  has  been 
optimized  for  use  on  a  massively  parallel  computer,  the  Thinking  Machine  CM-5,  that 
runs  at  MIT  as  part  of  an  ARPA  supercomputer  initiative,  project  SCOUT  (contract 
#MDA972-92-J-1032).  This  particular  machine  has  128  nodes,  each  with  a  single 
SPARC  microprocessor,  4  specialized  vector  units  and  32  MB  of  memory.  The  routine 
has  also  run  on  a  32  node  machine  at  NCAR. 

Massively  Parallel  computers  such  as  the  CM-5  consist  of  several  processing  units 
with  distributed  memory.  The  whole  system  is  interconnected  by  a  high-speed  com¬ 
munications  network.  Each  processing  unit  can  operate  independently,  but  total  pro¬ 
cessing  power  does  not  increase  linearly  with  the  number  of  processing  elements.  The 
speed  of  operations  depends  on  the  efficiency  with  which  information  is  supplied  to 
each  processor.  As  a  result,  different  types  of  communications  are  divided  by  efficiency 
into  a  hierarchy  operations  with  varying  desirability.  Even  with  a  high-speed  network, 
passing  information  between  distant  nodes  is  considerably  slower  than  accessing  near¬ 
est  neighbor  information,  which  in  turn  is  slower  than  reading  local  memory  units. 
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Thus  the  optimization  of  the  CM5  consists  of  coding  the  routine  such  that  it  minimizes 
all  inter-processor  communications. 

The  shallow  water  model  is  well-suited  to  parallel  computation  because  it  is  in¬ 
herently  local.  Each  processor  represents  a  grid  element,  and  the  time  stepping  of 
each  variable  depends  only  on  local  information  or  local  derivatives  which  are  nearest 
neighbor  operations.  Some  effort  is  required  to  insure  that  the  layout  of  the  many 
variables  are  aligned  properly  so  that  local  information  is  indeed  stored  locally.  When 
the  number  of  grid  elements  exceeds  the  number  of  processing  units,  the  CMS  compiler 
creates  virtual  units  that  automatically  optimize  local  information  transfers. 
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Chapter  3 


Linear  to  Nonlinear 
Cross-equatorial  Flow 


3.1  Introduction 

The  numerical  experiments  comprising  this  study  and  using  different  model  parame¬ 
ters  are  summarized  in  Table  3.1.  The  experiments  primarily  consider  upwelling  to  be 
uniform  following  Equation  2.14  with  a  few  runs  devoted  to  the  case  with  upwelling 
localized  at  high  latitudes  as  in  Equation  2.15.  The  linear  run  is  checked  against  the 
analytic  solution,  and  then  comparisons  are  made  with  more  strongly  forced  exper¬ 
iments.  Higher  Reynolds  numbers  axe  achieved  both  by  increasing  the  forcing  and 
keeping  dissipation  constant  and  by  maintaining  forcing  and  reducing  the  frictional 
boundary  layer  width.  An  upper  limit  on  the  forcing  is  that  which  produces  horizon¬ 
tal  velocities  approaching  the  gravity  wave  speed.  Model  runs  in  which  this  bound 
is  exceeded  fail  when  the  variable  h  vanishes  in  the  vicinity  of  the  equator  where  the 
rotational  constraint  is  weak,  and  thickness  advection  can  be  large  due  to  nonlinear 
steepening.  Some  comparison  is  made  between  the  runs  using  higher  forcing  values 
and  those  using  weaker  dissipation. 
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Run 

5o 

Ah 

Ho 

Si 

Rei 

Re2 

Comments 

U1 

.01 

.001 

2 

0.1 

0.022 

0.011 

5 

U2 

.04 

.001 

2 

0.1 

0.045 

0.089 

20 

U2.A 

.032 

.001 

1 

0.1 

0.056 

0.18 

32 

U3 

.1 

.001 

2 

0.1 

0.07 

0.35 

50 

U4 

.01 

.0002 

1 

0.058 

0.032 

0.16 

50 

U5.X 

.3 

.001 

2 

0.1 

0.12 

1.84 

150 

=  15 

U5.Y 

.3 

.001 

2 

0.1 

0.12 

1.84 

150 

Ly  =  30 

U6 

.2 

.0005 

2 

0.079 

0.1 

2.0 

200 

U7 

.2 

.00029 

2 

0.065 

0.1 

3.5 

350 

U8 

.1 

.0002 

1 

0.05 

0.1 

5.0 

500 

U9 

.05 

.00005 

1 

0.037 

0.07 

7.07 

1000 

=  0.0054 

UFS 

.2 

.0005 

1 

0.1 

0.1 

1.0 

200 

Free  Slip 

N2 

.1 

.00067 

1 

0.087 

0.07 

0.52 

75 

UBar 

lo-'* 

10~® 

1 

0.01 

0.014 

2.82 

100 

Table  3.1:  Table  of  numerical  runs  for  uniform  upwelling  configuration.  All  experiments  use  no 
slip  boundaries  except  for  that  noted.  Ho  refers  to  the  initial  undisturbed  depth  of  the  layer. 

3.2  Spin-up  of  the  Equatorial  Ocean 

Before  examining  the  steady-state  solution  it  is  useful  to  consider  the  spin-up  of  the 
system  from  a  state  of  rest.  The  evolution  of  the  flow  in  this  type  of  equatorial  ocean 
has  been  well-documented  by  Kawase  (1987).  That  work  used  linear  momentum 
dissipation  rather  than  the  Laplacian  form  used  here,  but  this  difference  is  minor  in 
the  limit  of  weak  flows.  Linear  wave  dynamics  describe  this  evolution  well. 

Illustrated  in  Figure  3.1  is  a  series  of  instantaneous  portraits  of  the  height  and 
velocity  fields  at  three  times  after  t  =  0,  when  the  forcing  is  initiated.  The  run 
is  frictional  and  forced  weakly,  with  5o  =  0.01,  Ah  =  10“^,  and  H  =  1.  Fluid 
entering  into  the  source  region  in  the  north-west  corner  of  the  basin  forms  a  small 
mound  which  begins  geostrophic  adjustment  within  a  rotation  period.  An  anti-cyclonic 
gyre  surrounding  the  source  forms  along  with  a  Kelvin  wave  that  travels  southward 
along  the  western  boundary.  In  its  wake  follows  the  fluid  entering  at  the  source. 
Upon  reaching  the  equator,  this  wave  couples  into  an  eastward  propagating  equatorial 
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Kelvin  wave.  Although  the  mass  source  is  located  in  the  northern  hemisphere,  the 
development  of  the  wave  field  is  nearly  symmetric  about  the  equator  soon  after  the 
signal  reaches  the  equator.  At  the  eastern  boundary,  the  disturbance  splits  into  two 
oppositely  directed  Kelvin  waves  travelling  along  that  boundary.  Finally,  anomalies 
at  the  eastern  boundary  accompanying  the  pa.ssage  of  each  Kelvin  wave  excite  Rossby 
waves  which  transmit  this  information  westward  across  the  interior.  By  establishing 
a  steady  resupply  of  fluid  from  the  western  boundary  current,  the  passage  of  the 
Rossby  wave  train  arrests  the  continued  drainage  of  the  interior  that  results  from  the 
up  welling  in  that  region.  Although  a  spectrum  of  Rossby  waves  is  generated  by  the 
Kelvin  wave  at  the  eastern  boundary,  in  practice  it  is  soon  after  the  passage  of  the 
fastest  mode  that  the  interior  reaches  an  approximately  steady  state.  At  the  western 
boundary,  a  steady  current  forms,  consistent  with  the  frictional  parameterization.  In 
the  present  study,  there  forms  a  southward  flowing  Munk  boundary  current  with  a 
weak  off-boundary  counterflow. 

The  final  state  of  the  system  includes  an  interior  flow  that  is  consistent  with  the 
linear  theory  of  Stommel  and  Arons  (1960a).  With  uniform  upwelling,  the  interior 
flow  is  weakly  poleward  and  eastward.  There  is  no  equatorial  signal.  Fluid  entering 
the  source-region  merely  follows  the  western  boundary  to  the  latitude  where  it  breaks 
from  the  coast  and  migrates  slowly  into  the  interior,  ultimately  to  be  upwelled  out  of 
the  basin. 

In  the  separate  case  of  a  localized,  southern  hemisphere  sink,  there  is  no  interior 
flow  in  the  final  state.  Two  intense  high-latitude  gyres  remain,  one  anti-cyclonic  for 
inflow,  and  one  cyclonic  for  outflow,  connected  by  the  boundary  current. 

The  development  of  this  flow  occurs  in  two  phases,  the  initial  Kelvin  wave  prop¬ 
agation  around  the  basin  boundary,  and  the  subsequent  Rossby  waves  which  return 
information  across  the  interior.  The  group  velocities  for  these  waves  indicate  the 
time-scale  for  spin-up. 
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1=10 


1=20 


t=30 


X 


Figure  3.1:  Time  sequence  of  height  contours  and  velocities  for  a  case  of  linear  spinup.  The  mean 
height  is  1,  with  contours  every  0.002.  The  initial  spinup  is  marked  by  the  passage  of  a  Kelvin  wave 
along  the  western  boundary,  the  equator,  and  the  eastern  boundary. 

The  Kelvin  waves,  which  move  at  the  gravity  wave  group  speed  of 


travel  about  the  basin  a  distance  L  —  Ly  +  L^-  Thus  the  time  scale  for  the  boundary 
field  is 


(3.2a) 

(3.2b) 


for  the  standard  run  of  this  study.  Figure  3.1c  shows  a  snapshot  at  this  time  for  a 
weakly  forced  run. 

The  Rossby  wave  time  scale  is  estimated  from  the  dispersion  relation 


(3.3) 
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The  fastest  westward-propagating  mode  is  that  for  which  horizontal  wavenumbers,  k 
and  /,  are  zero  (i.e.,  long  waves).  In  this  case,  the  group  velocity  is  given  by 


R 


(3.4a) 

(3.4b) 


in  the  scaling  of  this  model.  The  time  scale  for  Rossby  waves,  which  is  strongly  latitude 
dependent,  is  considerably  longer  than  for  Kelvin  waves.  With  a  maximum  latitude 
of  10  deformation  radii,  the  time  scale  for  adjustment  of  the  interior  is  approximately 


T 


R 


1000 


(3.5a) 


Thus  the  time-scale  for  Rossby  wave  propagation  is  significantly  longer  than  that 
for  the  Kelvin  waves  and  determines  the  required  time  for  the  linear  system  to  reach 
steady-state. 

One  measure  of  the  degree  to  which  the  system  has  reached  equilibrium  is  to 
compare  height  fields.  Figure  3.2  shows  the  difference  in  height  field  between  times 
1000  and  1200  for  the  linear  run,  Ul.  In  this  case  the  weak  up  welling  over  the  interior 
is  to*  =  5  X  10~®.  Thus  to* At  =  10“^.  The  variation  in  height  fields,  with  a  magnitude 
of  10“®,  has  some  structure  associated  with  the  slow  high  latitude  Rossby  waves  and 
very  weak  equatorial  waves,  but  the  system  is  very  near  to  a  steady  state. 

The  spin-up  of  a  strongly  forced  run  is  not  considerably  different  from  the  ca^e  de¬ 
scribed  above.  A  series  of  plots  presented  in  Figure  3.3  show  the  sequence.  The  initial 
descent  of  the  Kelvin  wave  along  the  western  boundary  is  similar,  and  it  couples  into 
an  equatorial  mode,  just  as  in  the  linear  case.  However,  the  group  velocity  comprises 
both  the  gravity  wave  speed  and  current  advection  [Philander  (1990)).  Furthermore, 
the  coupling  is  somewhat  more  complicated.  There  exists  an  inertial  overshoot  of  the 
Kelvin  wave  [Springer  and  Kawase  (1993)),  and  associated  with  that  cross-equatorial 
flow  is  an  instability.  As  a  result,  a  time  dependent  region  in  the  vicinity  of  the  western 


49 


boundary  and  equator  remains  behind,  even  after  the  Kelvin  wave  has  passed.  The 
subsequent  behavior  of  the  system  is  similar,  with  the  time-scale  to  reach  a  now  quasi¬ 
steady  state  largely  determined  by  that  of  the  Rossby  waves  in  the  interior.  Figure  3.4 
shows  the  difference  between  mean  height  fields  at  time  t  =  1100  and  t  =  3000  for  run 
U7.  Fields  are  averaged  over  100  time  units  to  remove  highest  frequency  fluctuations. 
In  this  case,  =  1,  and  the  average  height  differences  of  the  interior  are  a  fraction 

of  this  amount.  Only  near  the  dynamically  active  western  boundary  does  the  height 
difference  appear.  In  this  region,  low  frequency  fluctuations  remain  and  intense  fea¬ 
tures  shift  slightly  in  the  domain  over  long  time  periods.  In  addition,  this  region  is 
strongly  dominated  by  advective  processes  which  have  a  much  shorter  time-scale  than 
the  basin  waves.  As  in  the  linear  case  the  system  has  reached  a  quasi-steady  state. 
All  time-averaging  in  the  measurements  in  this  study  begin  at  time  t  =  1000  except 
run  U9  (which  begins  at  t  =  700)  in  which  the  zonal  basin  dimension  is  only  5  and 
therefore  the  spin-up  is  faster.  A  check  of  the  sensitivity  to  this  value  for  the  very 
nonlinear  run  (U8)  shows  very  minor  quantitative  differences  in  averaged  quantities 
that  begin  at  time  t  —  2000. 

The  time  dependence  is  transmitted  throughout  the  interior,  again  via  equatorial 
and  boundary  Kelvin  waves,  and  Rossby  waves,  but  the  height  variations  are  small 
compared  to  the  larger  anomaly  of  the  initial  pulse.  That  the  effect  of  nonlinearity 
is  not  more  strongly  felt  in  the  interior  is  a  testament  to  the  degree  to  which  linear 
physics  describe  fields  away  from  the  western  boundary  current  and  the  equator. 


3.3  Linear  Solution  and  Numerical  Experiment 

The  linear  solution,  which  can  be  determined  from  a  combination  of  Stommel  and 
Arons  (1960a)  theory  with  a  Munk  (1950)  boundary  layer,  provides  a  useful  baseline 
to  compare  the  more  complicated  nonlinear  flow.  This  section  reviews  these  two 
theories  using  the  no  slip  boundary  condition  only. 
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t=10  1=20  t=30 


3.3.  Time  sequence  of  height  contours  unci  velocities  for  a  case  of  a  nonlinear  spinup. 
The  undisturbed  height  is  1,  with  contours  every  0.1.  The  same  sequence  as  in  the  linear  case  is 
noted,  with  the  addition  of  remaining  variability  at  the  intersection  of  the  western  boundary  with 
the  equator.  The  Reynolds  number  for  this  run  is  1000. 
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Figure  3.4:  Difference  in  100  time  unit  averaged  height  fields  of  the  nonlinear  model  run  U7 
between  time  t  =  1100  and  t  =  3000  (xlO).  Some  structure  remains  at  the  western  boundary  where 
low  frequency  motions  shift  large  features,  but  the  system  has  reached  a  quasi-steady  state. 


Assuming  weak,  steady  currents  and  small  layer  height  variations  the  equatorial 
)0-plane  shallow-water  equations  reduce  to 


-Pyv 

fiyu 

Vu 


-9’ 


-gf 


d-q 

dx 

dq 

dy 


+  Ah 


d^v 

dx^ 


(3.5a) 

(3.5b) 

(3.5c) 


where  h  =  H  +  q,  and  q  H.  The  assumption  that  intense  boundary  currents  exist 
only  at  the  western  edge  of  the  basin  (which  results  from  the  existence  of  nonzero 
P)  has  been  made  implicitly  in  the  one  component  of  viscosity  that  remains  in  the 
equations.  The  meridional  velocity  is  in  geostrophic  balance. 

The  linear  vorticity  equation,  which  follows  directly  from  Equations  3.5a-3.5c,  ex¬ 
presses  a  balance  of  planetary  advection  of  vorticity  with  stretching  and  vorticity 
diffusion: 

No  slip  and  no  flux  boundary  conditions  apply  along  the  basin  boundary: 


u  =  0;u  =  0 


(3.7) 


It  is  convenient  to  separate  variables  into  interior  and  boundary  layer  components. 


u  =  ui  +  ub]  v  =  vi  +  vb]  q=qi  +  qB-  (3.8) 

Taking  the  Munk  boundary  layer  parameter,  Sm  L^,  the  interior  fields  are  deter¬ 
mined  to  be 


VI 


ui 


hi 


w^. 


=  y 


=  / 


=  H 


H 

Xe 


2w 


■*  y 

H  ^  H  dy 


-j: 


Py^w* 

g'H 


dx 


dx 


(3.9a) 

(3.9b) 

(3.9c) 
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The  boundary  layer  structure  results  from  a  balance  of  planetary  vorticity  advec- 
tion  and  vorticity  diffusion.  This  balance  is  expressed  in  terms  of  the  layer  height 


variable,  tjb, 


Thus 


dvB  ^  p  d'^VB 
dx  ^  dx^ 


(3.10) 


TjB  =  Cicos  +  C2sin  (3  H) 

Variables,  Ci,  are  functions  of  latitude,  y.  That  the  boundary  layer  variables  must 
vanish  as  ^  becomes  0(Lx)  requires  C3  =  Ca  =  0.  The  no  slip  condition  sets  to 
leading  order 


(3.12) 


The  no  flow  condition  determines  C2  by  integrating  the  height  equation.  Equa¬ 
tion  3.5c,  over  a  sub-area.  A,  of  the  basin,  bounded  to  the  north,  east,  and  west  by 
basin  boundaries,  and  to  the  south  by  a  line  of  latitude,  as  in  Figure  3.5.  Given  no 
flux  and  no  slip  boundary  conditions,  this  integrated  equation  reduces  to 


(3.13) 


and  simply  restates  the  conservation  of  mass  within  the  domain,  A. 

Inserting  the  geostrophic  relation  for  v  and  separating  interior  and  boundary  layer 
components  gives 


j3y  f^E  r  rVN 

Cl  =  VB{0,y)  = - —  /  /  w^dy  -  yw^  . 

Q  II  Jo 


(3.14) 


Horizontal  velocities  for  the  boundary  layer,  ub  and  vb,  follow  directly  from  the  gov¬ 


erning  equations. 


g'  dr}B 

fiy  dx 

d^VB 

Py  dy  dx^ 


(3.15a) 

(3.15b) 
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Figure  3.5:  Subdomain  A  of  full  basin  bounded  by  a  border  T  comprising  the  eastern,  western, 
and  northern  boundaries,  as  well  as  a  line  of  latitude  at  y  =  j/q. 


Although  the  fields  depend  strongly  on  the  distribution  of  the  forcing,  tn*,  some 
features  remain  robust.  The  interior  is  locally  forced,  with  the  meridional  velocity 
linearly  proportional  to  the  local  upwelling.  The  zonal  interior  velocity  field  acts  to 
distribute  mass  longitudinally  to  accommodate  the  upwelling  and  meridional  velocity. 
The  interior  height  field  simply  reflects  a  geostrophic  balance  with  the  velocity  fields. 
It  is  helpful  to  consider  the  simplest  of  cases,  in  which  the  upwelling  is  constant, 
=  Wo,  and  there  exists  a  source  of  fluid  having  a  mass  flux,  5o,  entering  through 
the  northern  boundary.  The  dynamics  associated  with  the  mass  flux  will  be  ignored, 
except  in  that  it  provides  the  necessary  mass  to  the  system  to  be  expelled  through  the 
layer  interface.  Under  this  assumption. 


VI  = 

ui  = 
hi  — 


Wo 
^  H 
2Wol. 

H 

g'H 


(3.16a) 

(3.16b) 

(3.16c) 
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where,  =  xe  —  x,  and 


C'l  =  (2y  -  vn)  -  5o] .  (3.17) 

g'h 

As  before  L*  is  the  full  zonal  scale  of  the  basin.  This  case  is  particularly  simple, 
and  both  the  height  field  and  the  horizontal  velocities  are  illustrated  in  Figure  3.6. 
Parameters  for  this  plot  are  =  10,  ytf  =  10,  Wq  =  5  x  10~®,  =  —WoL^Ly, 

g'  =  =  I  and  the  undisturbed  height  is  2H  The  Munk  boundary  layer  width, 

6m  =  0.1.  Height  contours  are  every  0.25H  except  for  1.9  <  h  <  2.1,  where  the  interval 
is  0.025H.  This  figure  shows  the  classic  Stommel  and  Arons  interior  solution  with  a 
Munk  boundary  layer.  The  interior  has  weak  poleward  and  eastward  flow  in  both 
hemispheres.  The  velocity  vectors  are  naturally  small  in  regions  where  the  flow  is  weak 
and  difficult  to  distinguish  in  the  figure.  In  these  regions,  the  sense  of  the  flow  must 
be  inferred  from  the  height  contours.  There  is  an  intense  western  boundary  current 
that  brings  fluid  south  from  its  entrance  into  the  basin  at  the  northern  boundary. 
There  appears  to  be  considerable  flux  through  the  northern  and  southern  boundaries, 
both  in  the  interior  and  the  western  boundary.  It  is  possible  to  extend  the  above 
calculations  to  include  northern  and  southern  boundary  currents,  but  these  features 
are  not  particularly  illuminating,  except  that  they  transfer  mass  back  from  the  interior 
of  the  basin  to  the  western  boundaries. 

One  characteristic  feature  of  the  Stommel  and  Arons  (1960a)  theory  is  that  the 
transport  of  the  boundary  current  vanishes  north  of  the  southern  boundary  of  the 
basin.  In  the  figure,  it  is  apparent  that  this  occurs  at  y  =  —5,  and  this  latitude  can  be 
calculated  from  the  mass  balance  of  the  interior  flow.  Consider  a  basin  whose  northern 
boundary  coincides  with  the  equator  and  whose  source  is  located  at  this  northern  edge. 
The  mass  balance  at  any  line  of  latitude  yields 

f^E  fy  f^B  to 

/  {vi  +  VB)dx-  /  —dxdy.  (3.18) 

Jo  J J 0 

Inserting  Equation  3.9a  and  setting  the  boundary  layer  transport  to  zero  determines 
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Figure  3.6:  Steady  solution  to  Equation  3.6  with  uniform  upwelling  and  a  Munk  boundary  layer 
width  of  0.1. 


the  critical  latitude.  In  the  case  when  upwelling  is  uniform, 

1 

Vc  =  -ys- 


When  upwelling  is  localized  as  a  step  function  south  of  latitude,  t/i, 

Vc  =  ^(^5  -yi). 


(3.19) 

(3.20) 


In  the  current  example,  ys  =  —10  and  =  0  which  gives  the  appropriate  value 
for  yc-  This  location  divides  the  boundary  current  into  two  halves,  one  driven  by  the 
existence  of  the  northern  source  feeding  the  basin  upwelling,  and  a  second  current, 
driven  by  the  local  upwelling  of  the  southern  hemisphere,  independent  of  the  source  in 
the  north.  Fluid  destined  to  reach  the  southern  portion  of  the  basin  near  the  western 
boundary  has  an  extensive  journey  from  the  western  boundary  in  the  north  through 
the  interior  in  the  south  and  back  to  the  western  boundary  before  exiting  the  basin. 

The  time-dependent  numerical  model  has  been  run  to  test  the  code  in  this  linear 
limit.  Shown  in  Figure  3.7  are  the  height  contours  and  horizontal  velocities  of  a 


58 


Figure  3.7:  Steady  state  achieved  by  the  numerical  model,  using  uniform  upwelling,  except  for 
the  source  region  in  the  north  west  corner.  Compare  to  Figure  3.6.  Note  that  the  Reynolds  number 
for  this  run  is  Jie  =  5. 

numerical  model  run,  Ul,  that  corresponds  directly  to  that  of  the  analytic  calculation 
above.  In  the  north-western  corner  of  the  basin  is  a  region  of  large  velocities  and 
height  gradients.  This  is  the  source  region  for  fluid  in  the  numerical  model,  and 
its  dynamics  are  not  included  in  the  theory  presented  above.  Also  visible  are  the 
northern  and  southern  boundary  currents  that  return  mass  to  the  western  boundary. 
The  contours  in  both  flgures  are  identical  and  can  be  compared  directly.  The  numerical 
and  ajialytically  calculated  fields  match  well,  far  from  the  region  of  strong  forcing  and 
the  northern  and  southern  boundaries. 

As  expected,  the  system  divides  into  multiple  gyres,  one  very  small  and  associated 
with  the  localized  source  of  fluid  to  the  system  and  two  hemisphere-scale  gyres  divided 
by  the  equator.  In  this  linear  limit,  the  potential  vorticity  contours,  not  shown,  simply 
follow  lines  of  latitude. 
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3.4  Transition  to  Instability 

As  the  Reynolds  number  increases,  the  steady  linear  solution  becomes  deformed  by 
the  influence  of  nonlinearity  in  the  dynamics.  The  interior  remains  largely  linear 
in  its  balance,  but  inertial  effects  become  important  in  the  boundary  layer.  As  the 
meridional  velocity  of  the  boundary  current  increases,  the  relative  vorticity  in  that 
region  becomes  non-negligible.  The  qualitative  distinction  from  the  linear  solution 
is  that  the  potential  vorticity  contours  now  detour  from  /-contours  near  the  western 
boundary.  When  the  boundary  condition  is  no  slip,  the  PV  contours  shift  northward 
just  at  the  boundary.  When  the  condition  is  free  slip,  the  contours  shift  only  slightly  at 
the  boundary,  the  result  of  small  variations  in  the  layer  height,  but  more  significantly, 
they  extend  south  of  the  lines  of  latitude  just  off  the  boundary. 

At  a  critical  Reynolds  number,  20  <  Rec  <  50,  the  system  becomes  unstable,  and 
a  steady  state  solution  is  no  longer  achieved  by  the  time- dependent  model.  However, 
a  quasi-steady  state  is  reached.  Instantaneous  pictures  of  the  height  and  velocity 
fields  show  the  development  of  periodic  eddies  just  north  of  the  equator  in  the  western 
boundary  current.  An  example  from  Run  U3  with  Re  =  50  is  presented  in  Figure  3.8. 
Only  50  vectors,  linearly  interpolated  from  the  stretched  onto  the  regular  grid  shown, 
span  the  basin  in  the  figure.  However,  approximately  50  grid  points  in  the  model 
successfully  resolve  eddies  in  this  region. 

Formed  north  of  the  equator,  these  eddies  are  always  cyclonic,  reflecting  the  high 
potential  vorticity  of  their  source  fluid  relative  to  the  local  planetary  vorticity.  The 
eddies  migrate  coherently  southward  along  the  boundary,  dissipating  vorticity  along 
their  way  by  creating  anomalously  large  velocity  gradients  in  the  vicinity  of  the  bound¬ 
ary.  As  the  eddies  propagate  across  the  equator  their  vorticity  decreases  until  they 
join  smoothly  with  the  steady  southern  hemisphere  flow  to  the  south. 

This  state  of  the  flow  is  time-dependent,  but  still  only  weakly  nonlinear.  The 
system  has  a  periodicity  associated  with  these  features,  and  their  geographic  extent 
is  extremely  limited,  one  eddy  length  scale  east  into  the  basin,  and  a  few  deformation 
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Figure  3.8:  Instantaneous  height  contours  and  velocity  vectors  for  run,  U3.  Re  =  50. 
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radii  north  and  south  of  the  equator.  The  eddy  motion  is  very  regular  and  time-series 
analyses  (not  shown)  of  the  potential  vorticity  at  various  locations  in  the  basin  indicate 
periodic  structure.  The  time-scale  associated  with  the  eddies  is  0(20  which 

corresponds  to  0(40  days)  in  dimensional  units.  This  time-scale  can  be  adjusted  by 
tuning  parameters  in  the  model,  but  realistic  scalings  for  the  transport  and  layer 
depth  yield  time-scales  that  compare  sensibly  with  observations  of  variability  as  in 
Richardson  and  Schmitz  (1993)  or  Johns  et  al.  (1993).  A  more  detailed  discussion  of 
the  eddy  generation  process  is  found  in  Chapter  5. 

Additional  information  can  be  gained  by  considering  time-averaged  fields.  Fig¬ 
ure  3.9  shows  the  height  contours  and  velocity  field  for  the  basin  averaged  over  2000 
time  units  of  the  model  (0(100)  eddy  time-scales)  after  spin-up  at  t  =  1000 
In  the  mean,  these  fields  look  quite  similar  to  the  linear  solution.  The  height  contours 
are  deformed  somewhat  and  the  potential  vorticity  contours  are  more  exaggerated  at 
the  western  boundary,  but  overall  there  is  little  qualitative  difference. 

3.5  Strongly  Nonlinear  Flow 

At  still  larger  Reynolds  number,  the  influence  of  nonlinearity  becomes  ever  stronger. 
The  eddies  grow  more  intense  and  strongly  deform  mean  potential  vorticity  contours. 
The  height  contours  extend  farther  to  the  south.  The  geographic  extent  of  the  vari¬ 
ability  extends  over  a  broader  portion  of  the  basin,  both  meridionally  and  zonally 
along  the  equator. 

Figure  3.10  shows  the  instantaneous  height  and  velocity  fields  for  a  run  at  Re  =  500 
at  time,  t  ~  2900.  The  eddy  field  is  extremely  disordered,  populating  a  broad  swath 
of  the  western  part  of  the  basin.  Eddies  extend  south  of  the  critical  latitude,  y^, 
where  boundary  currents  meet.  There  is  a  clear  indication  of  wave  propagation  out 
along  the  equator  and  eastern  boundary,  and  some  variability  is  visible  in  the  interior 
as  deformations  of  the  height  contours  in  that  region.  The  existence  of  such  strong 
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flows  near  the  equator  suggests  a  possible  mechanism  to  generate  the  tracer  tongue 
of  chloroflouromethane  observed  by  Weiss  et  al.  (1985)  and  shown  in  Figure  1.4.  A 
calculation  of  the  time-mean  velocity  (shown  below)  reveals  a  meandering  current 
near  the  equator  and  the  western  boundary  that  weakens  quickly  eastward.  More 
likely  responsible  for  the  tongue  is  the  Stokes  drift  associated  with  the  time-varying, 
equatorial  waves  visible  in  Figure  3.10.  This  possibility  was  first  suggested  by  Kawase 
and  Sarmiento  (1986)  and  seems  plausible  given  the  strong  generation  of  equatorial 
waves  in  this  model. 

The  instantaneous  potential  vorticity  field  provides  a  similarly  complicated  illus¬ 
tration  of  the  anomalies  in  this  highly  nonlinear  run.  This  field  is  presented  in  Fig- 
which  is  shaded  to  aid  the  distinction  between  positive  (light)  and  negative 
(dark)  values.  Direct  correlation  between  strong  features  in  the  two  figures  is  evident 
in  isolated  regions  of  anomalous  potential  vorticity.  These  regions  have  various  di¬ 
mensions,  though  a  few  structures  near  the  western  boundary  have  a  common  scale 
of  about  1  deformation  radius.  Eddies  of  both  sign  are  now  created  in  the  formation 
process  and  are  found  in  either  hemisphere,  though  generally  not  far  from  the  equator. 
In  addition,  there  exist  large  gradients  in  the  potential  vorticity  field,  particularly  near 
the  coherent  eddy  of  the  northern  hemisphere  and  near  the  western  boundary.  Indeed 
it  is  not  difficult  to  imagine  that  these  large  gradients  are  crucial  in  the  process  of 
potential  vorticity  transformation  as  will  be  seen  in  the  next  chapter. 

Also  visible  in  Figure  3.11  is  the  potential  vorticity  of  the  source  fluid.  Since 
this  quantity  depends  on  the  dynamical  behavior  of  the  source  region  it  is  not  easily 
related  to  the  input  parameters  of  the  system  and  must  be  diagnosed  from  the  results. 
The  source  fluid  in  this  experiment  has  a  value  of  potential  vorticity  between  4  and 
6  {(3LdH)~^,  though  the  fluid  in  the  narrow  western  boundary  current  has  a  range  of 
values.  In  the  northern  hemisphere,  negative  potential  vorticity  is  found  immediately 
adjacent  to  the  boundary  because  of  the  extreme  shear,  and  the  zonal  gradient  is 
very  large,  with  potential  vorticities  of  magnitude  6  just  a  fraction  of  a 


I  I  I  >  I  -  r - 1  r: - - - - — c — i 

0123456789  10 


X 

Figure  3.10:  Instantaneous  fields  for  very  nonlinear  run  (Re  =  500)  at  time  t  =  2900.  Cyclonic 
and  anti-cyclonic  eddies  of  various  scales  populate  the  basin. 
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deformation  radius  away. 

As  with  the  run  at  intermediate  Reynolds  number,  the  time-averaged  fields  provide 
some  order  to  this  complicated  tangle.  Shown  in  Figure  3.12,  the  time-averaged  height 
field  reveals  a  mean  cross  equatorial  flow,  just  as  in  the  low  Reynolds  number  runs, 
with  a  very  intense  western  boundary  current.  However,  distinguishing  this  run  is  the 
large  feature  that  resides  just  north  of  the  equator.  This  large  cyclonic  eddy  enhances 
the  transport  of  the  western  boundary  current  and  redirects  a  portion  of  the  fluid  out 
into  the  interior  and  back  to  the  northern  hemisphere  where  it  migrates  its  way  into 
the  interior.  The  magnitude  of  this  feature  is  directly  linked  to  the  substantial  inertial 
overshoot  of  fluid  in  the  western  boundary  current. 

To  further  understand  the  behavior  of  the  mean  field,  it  is  helpful  to  examine  the 
meridional  mass  transport  within  the  basin,  integrated  from  the  western  boundary; 

Tm  =  hvdx.  (3.21) 

This  quantity  would  be  equivalent  to  the  mass  transport  stream  function  were  the 
system  horizontally  non-divergent.  Since  the  experiment  includes  weak  upwelling 
throughout  the  domain,  it  is  necessary  to  consider  this  quantity  the  meridional  trans¬ 
port  only.  In  regions  of  intense  western  boundary  flow,  streamlines  and  meridional 
transport  contours  must  be  nearly  coincident,  but  not  in  regions  of  the  interior. 

In  Figures  3.13  and  3.14  the  function  Tm  is  contoured  for  both  low  and  high 
Reynolds  number  flows.  The  volume  of  fluid  passing  through  the  basin  is  Q.IHcLd 
in  both  cases,  and  half  of  the  volume  is  upwelled  in  either  hemisphere.  Observe  that 
in  the  more  frictional  case.  Figure  3.13,  the  —0.075  contour  extends  quite  far  south, 
nearly  to  the  equator.  Almost  three  quarters  of  the  fluid  entering  the  basin  follows 
the  western  boundary  to  the  equator.  The  contrasting,  more  nonlinear  case  shown 
3.14  reveals  that  fully  all  of  the  water  that  eventually  will  upwell  in  the 
northern  hemisphere  of  the  model  extends  down  to  the  equator. 

It  is  important  to  distinguish  this  eddy-like  feature  of  the  mean  field  from  a  station¬ 
ary,  recirculating  eddy.  The  appearance  of  closed  contours  in  the  mean  suggests  that 
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Figure  3.11:  Instantaneous  PV  field  for  very  nonlinear  run  {Re  =  500)  at  time  t  =  2900.  Positive 
potential  vorticity  has  the  lighter  shades  and  negative  darker.  Pure  white  corresponds  to  negative 
values  less  than  12. 
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Figure  3.13:  Time-averaged  meridional  transport  for  the  Re  =  50  run.  The  volume  flux  through 
the  source  and  sink  is  (SAHcLd, 


Figure  3.14:  Time-averaged  meridional  transport  for  the  Re  =  500  run,  Sq  =  QAEcLd- 
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water  is  trapped  in  this  region  just  north  of  the  equator.  To  the  contrary,  the  coherent 
structures  in  the  instantaneous  snapshots  grow,  migrate  southward,  and  decay  con¬ 
tinually.  I  believe  the  appearance  of  the  mean  eddy  is  an  indication  that  the  location 
simply  experiences  protracted  eddy  growth  or  sustenance  before  it  ejects  anomalies 
southward  across  the  equator. 

One  indication  that  there  is  no  trapped  fluid  within  this  recirculation  is  found  in 
the  time-averaged  potential  vorticity  field  shown  in  Figure  3.15.  There  is  no  indication 
of  isolated  patches  of  anomalous  potential  vorticity  in  the  region  of  the  mean  eddy, 
and  it  is  unlikely  that  potential  vorticity  would  be  shifting  continually  back  and  forth 
between  portions  of  the  recirculation.  Thus  the  mean-eddy  does  not  reflect  a  local 
stationary  recirculation  in  the  sense  of  trapped  fluid. 

Also  visible  in  Figure  3.15  is  the  strong  gradient  of  potential  vorticity  at  the  western 
boundary.  Indeed,  contours  in  the  mean  remain  considerably  close  to  /-contours, 
except  in  the  immediate  vicinity  of  the  western  boundary  where  height  variation  and 
particularly  velocity  shear  are  extreme.  With  no  slip  boundary  conditions,  the  zero 
potential  vorticity  contour  extends  several  deformation  radii  up  the  western  boundary. 
As  a  reality  check,  it  is  helpful  to  reconsider  the  observations  of  mean  velocity  in  the 
Atlantic  ocean  shown  in  Figure  1.3.  Estimating  the  relative  vorticity  in  the  strong 
shear  zone  suggests  that 

V  .5ms~^ 

Sx  50km 

-  10-^5-^ 

Values  of  the  planetary  vorticity  scale  with 

PLd  =  2  X  10““m“^s~^  ■  250ifcm 
^  5  X  10-®s-\ 

Thus  the  prospect  of  negative  values  of  mean  planetary  vorticity  creeping  up  the 
western  boundary  is  not  inconsistent  with  data. 
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Figure  3.15:  Time-mean  potential  vorticity  for  Run  U8  {Re  =  500). 
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3.6  Summary 


The  numerical  model  successfully  transports  fluid  across  the  equator  for  a  range  of 
Reynolds  numbers.  At  low  Re,  cross-equatorial  flow  is  laminar,  and  the  flow  is  well 
predicted  by  linear  theory.  At  a  critical  Reynolds  number,  the  boundary  layer  becomes 
unstable,  and  the  model  creates  cyclonic  eddies  which  propagate  periodically  across 
the  equator.  Very  high  Reynolds  number  flow  is  characterized  by  considerably  more 
variability,  with  cyclonic  and  anti-cyclonic  features  covering  an  extensive  portion  of 
the  tropical  ocean,  particularly  near  the  western  boundary  and  equator.  Regardless 
of  strength  of  the  nonlinearity,  the  system  adjusts  so  as  to  allow  cross-equatorial 
flow.  The  subject  of  the  next  chapter  is  the  characterization  of  the  potential  vorticity 
transformation  in  this  flow. 
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Chapter  4 


Vorticity  Analysis  and  Sensitivity 


This  chapter  returns  to  the  dynamical  problem  of  vorticity  transformation  in  inertio- 
frictional  cross-equatorial  flow  and  quantitatively  analyzes  the  vorticity  fluxes  in  the 
set  of  numerical  experiments  described  in  Chapter  3.  The  first  section  defines  the  J- 
vector  as  the  primary  tool  of  the  analysis  and  examines  the  fluxes  for  flows  of  different 
Reynolds  number.  The  dynamical  balances  of  the  time-mean  western  boundary  layer 
is  subsequently  decomposed. 

At  the  end  of  the  chapter  are  sections  that  test  the  sensitivity  of  the  vorticity 
fluxes  to  different  model  configurations.  In  particular,  the  effect  of  free-slip  boundary 
conditions  and  non-uniform  upwelling  distributions  are  considered. 


4.1  Time- averaged  Vorticity  Analysis 

4.1.1  J- vectors 

Definitions 

By  averaging  Equation  2.16  over  time-scales  that  are  long  compared  to  the  time-scales 
for  fluctuations  in  the  system,  one  obtains  the  steady-state  vorticity  balance.  Letting 
overbars  with  a  t  superscript  denote  time  averaging,  the  steady-state  vorticity  equation 
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becomes 


u  •  vc‘  +  tj*  +  (y  +  C)V  •  u‘  =  ~AhV^(()\  (4.1) 

expressing  the  balance  of  of  four  components,  the  nonlinear  advection  of  relative  vor- 
ticity,  advection  of  planetary  vorticity,  stretching  (which  includes  a  nonlinear  compo¬ 
nent),  and  friction. 

Equation  4.1  can  be  further  simplified  by  expressing  it  in  flux  form.  In  terms  of 
the  absolute  vorticity,  =  y  -f  Equation  4.1  is  restated,  dropping  the  superscripts. 


V  ■  (uCx)  =  AhV2(C). 

(4.2) 

Rearranging  terms. 

o' 

II 

> 

(4.3) 

where 

j  —  —  A;fV((^) 

(4.4a) 

=  ua  +  u'G-AHV(0 

(4.4b) 

Here,  the  primed  variables  represent  deviations  from  their  time-mean  quantities. 
Physically,  J-vectors  represent  horizontal  vorticity  fluxes  within  the  fluid  layer.  There 
is  zero  cross-isopycnal  flux,  and  the  field  is  horizontally  non-divergent.  Although 
the  significance  of  the  vorticity  equation  has  an  extensive  history  dating  back  to  the 
last  century,  the  J-vector  formulation  for  three  dimensional  systems  is  more  recent 
{Haynes  and  McIntyre  (1987)).  In  that  work,  J-vectors  depict  the  vorticity  fluxes 
along  isentropic  sheets,  whose  relationship  to  a  single  layer  model  follows  directly 
from  an  isentropic  coordinate  representation.  Indeed  the  shallow-water  model  can  be 
considered  to  represent  a  single  isentropic  layer  of  the  ocean. 

It  IS  convenient  to  divide  the  vorticity  flux  into  three  components,  one  due  to  mean 
advection  of  mean  absolute  vorticity,  another  to  the  eddy  flux  of  absolute  vorticity, 
and  the  third  resulting  from  the  flux  of  vorticity  due  to  frictional  forces  acting  in  the 
system: 

J  ^mn  T  Jed  4"  J/rtc  •  (4.5) 
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The  simplicity  of  the  J-vector  notation  is  now  evident.  By  considering  the  mass- 
weighted  potential  vorticity  rather  than  the  full  potential  vorticity,  the  forcing  in  the 
height  equation  is  incorporated  implicitly  in  the  analysis.  As  a  result,  only  three 
terms  are  responsible  for  a  dynamical  balance,  as  opposed  to  five  in  the  full  potential 
vorticity  equation  when  similarly  decomposed.  Of  course,  this  formulation  has  its 
cost:  the  loss  of  any  dynamical  information  associated  with  the  forcing  in  the  height 
equation. 

Since  J  is  non-divergent,  the  flux  entering  any  region  must  balance  with  flux  leav¬ 
ing.  Integrating  Equation  4.3  over  any  area,  A,  bounded  by  a  border,  F,  yields 

jfj-n  =  0.  (4.6) 

The  numerical  approximation  to  the  J- vectors  specific  to  the  discretization  is  de¬ 
scribed  in  Appendix  A.4. 


4.1.2  Analysis  and  Numerical  Results 


Linear  Analysis 

Before  proceeding  to  the  numerical  J-vector  plots,  it  is  instructive  to  consider  re¬ 
sults  for  linear  flow  first.  When  the  dynamics  are  linear,  the  J-vectors  reduce  to  the 
following. 


Jed  =  0 


^  fric 


yu 


=  -Ah 


'  d^v  d^v 
dxdy  ) 


{4.7a) 

(4.7b) 

(4.7c) 


Embedded  in  the  above  approximations,  are  the  assumptions  that  the  flow  is  steady 
and  that  relative  vorticity  is  small  compared  to  planetary  vorticity  and  the  boundary 
layer  approximation  ^  ^ 

It  is  important  to  note  that  the  advective  meridional  vorticity  flux  at  any  location 
is  only  a  function  of  the  meridional  velocity  and  the  latitude:  =  vy.  As  a  result. 
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the  meridional  flux  of  J,  while  non-zero  at  all  mid-latitudes  with  non-zero  v,  must 
vanish  at  y  =  0.  At  the  equator  in  the  linear  model,  there  can  be  no  advective  flux  of 
vorticity  at  any  longitude,  regardless  of  the  meridional  velocity.  The  meridional  flux 
due  to  friction, 

'^frlc  =  -^nCy,  (4.8) 

is  not  necessarily  zero  at  the  equator.  However,  integrating  this  quantity  from  the 
western  boundary  gives 


-  Ah  f  v^ydx  =  ~Ah  (vy{x)  -  ny(0))  .  (4.9) 

The  rightmost  term  on  the  right-hand  side  is  zero  with  no  slip  boundary  conditions, 
and  the  remaining  contribution  is  small  if  the  integration  is  carried  into  the  interior.  If 
the  boundary  layer  structure  is  independent  of  y  and  the  interior  fields  are  symmetric 
about  the  equator,  then  the  meridional  flux  of  vorticity  due  to  friction  is  exactly  zero. 
Relaxing  these  constraints,  but  integrating  across  the  domain,  x  —  xh^  shows  that  the 
net  integrated  frictional  flux  of  vorticity  at  any  latitude  is  in  fact  zero.  The  remaining 
discussion  ignores  any  contribution  to  the  meridional  flux  due  to  friction. 

Recall  the  Stommel  and  Arons  flow  with  simple  Munk  boundary  layer  dynamics 
described  in  section  3.3.  The  meridional  flow  in  the  interior  vanishes  at  the  equator, 
but  the  boundary  layer  flow  remains.  It  carries  the  necessary  mass  across  the  equator 
to  upwell  in  the  southern  hemisphere. 

Equation  4.6  requires  that  the  net  integrated  flux  of  vorticity  into  a  sub-domain  of 
the  basin  that  is  bounded  by  the  equator,  the  western  and  eastern  boundaries  and  a 
line  of  latitude,  yo,  north  of  the  equator,  must  be  zero.  The  vorticity  flux  across  each 
portion  of  the  boundary  can  be  considered  independently.  From  the  north. 


'  yovdx 
0 

y^Northj 


(4.10a) 

(4.10b) 


where  Vnotik  is  the  total  mass  transport  per  unit  depth  through  the  northern  boundary 
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of  the  sub-domain.  Since  the  net  transport  across  this  latitude  is  southward,  there  is 
a  net  vorticity  flux  into  the  region  across  this  boundary. 

Since  the  equator  is  a  boundary  of  vorticity  flux,  there  is  no  contribution  to  sub- 
domain  through  the  southern  boundary.  In  addition  there  is  no  flow  through  the 
eastern  boundary  and  no  eastern  boundary  current  which  might  result  in  a  frictional 
flux  there.  As  a  result,  the  vorticity  flux  through  the  eastern  boundary  is  zero  also. 
The  only  sink  for  iNorth  is  the  western  boundary,  and  it  results  from  the  strong  shear 
in  the  boundary  current: 

= -A„  £  ^^dy.  (4.11) 

Thus,  vorticity  enters  the  region  through  the  northern  boundary  and  leaves  through 
the  western  boundary,  and  this  pathway  resides  entirely  in  the  northern  hemisphere. 

It  is  possible  to  further  dissect  the  structure  of  the  vorticity  fluxes  by  considering 
the  structure  of  the  velocity  field.  Although  the  net  transport  at  middle  latitudes 
is  negative,  in  the  interior  there  exists  a  slow  poleward  velocity,  and  therefore  an 
accompanying  weak  poleward  vorticity  flux  (i.e.,  out  of  the  domain).  This  poleward 
flux  is  compensated,  and.  sometimes  overcompensated,  by  the  equatorward  flux  in 
the  intense  western  boundary  current.  The  domain  bounded  by  the  equator,  the 
meridional  boundaries  and  the  line  of  latitude,  y,  can  be  dissected  by  a  line,  xo{y), 
such  that 

nvdx  =  0.  (4.12) 

Jxq 

All  J  entering  the  region  ea.st  of  Xq  leaves  again  through  the  northern  boundary  of  the 
subdomain.  All  J  entering  west  of  xq  is  expelled  through  the  western  boundary. 

A  schematic  diagram  of  the  vorticity  fluxes  under  linear  dynamics  is  shown  in 
Figure  4.1.  The  diagram  is  somewhat  exaggerated  in  the  zonal  scale  in  order  to 
broaden  the  western  boundary  current  with  respect  to  the  basin  dimension.  The 
strong  southward  flux  at  any  latitude,  y  =  t/o,  is  marked  near  the  western  boundary, 
as  is  the  bounding  vorticity  streamline,  x  —  xo{y).  The  vorticity  flux  west  of  xq  must 
exit  the  basin  at  the  western  boundary  whereas  the  vorticity  flux  immediately  to  its 
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east  is  exported  by  the  northward  Sverdrup  flow  of  the  interior. 

For  completeness,  it  is  useful  to  understand  J-vectors  south  of  the  equator.  The 
same  arguments  above  apply  to  a  sub-domain  that  is  bounded  to  the  north  by  the 
equator  and  to  the  south  by  a  line  of  latitude.  With  y  <  0,  the  vorticity  fluxes  due 
to  the  mean  field  are  reversed,  though  the  frictional  flux  is  not.  Thus  south  of  the 
equator,  the  J  vector  diagram  is  inverted  from  its  northern-hemisphere  counterpart. 

The  J- vector  quiver  plot  for  a  numerical  experiment  with  i2e  =  5  (Ul)  is  presented 
for  comparison  in  Figure  4.2.  Only  a  western  subdomain,  with  x  =  [0,2]  and  y  = 
[—5,5],  of  the  basin  is  illustrated  to  emphasize  the  fluxes  associated  with  the  western 
boundary  current.  An  equatorword  vorticity  flux  is  visible  near  the  western  boundary 
in  both  hemispheres  veering  westward  and  out  the  western  boundary  by  frictional 
processes.  In  addition,  a  portion  of  the  vectors  in  the  eastern  part  of  the  boundary 
current  turn  eastward  toward  the  interior  where  a  weak  eastward  and  northward  flux 
of  vorticity  is  associated  with  the  interior  flow.  These  vectors  are  difficult  to  see  in 
the  diagram  due  to  their  small  magnitude. 

An  alternative  presentation  reveals  more  clearly  the  vorticity  fluxes  in  the  basin. 
The  horizontal  non-divergence  of  the  vector  field  permits  the  construction  of  a  po¬ 
tential  function,  similar  to  the  streamfunction  for  a  two-dimensional  velocity  field 
such  that 

J  =  k  X  (4.13) 

Contours  of  this  field  are  streamlines  of  the  vorticity  flux;  lines  of  constant  'J  are 
tangent  to  the  J-vectors.  An  example  for  the  linear  run  U1  is  presented  in  Figure  4.3. 
Although  the  vorticity  flux  itself  on  any  given  line  is  not  constant,  contours  of  equal 
interval  in  ^  demarcate  tubes  of  constant  vorticity  flux.  Regions  of  large  gradients 
in  ^  (i.e.,  with  compressed  contours)  represent  regions  of  intense  vorticity  transport, 
and  conversely,  weak  gradients  correspond  to  weak  flux  of  J. 

In  Figure  4.3,  the  strong  equatorward  flux  of  vorticity  is  visible  near  the  western 
boundary.  The  trend  of  some  vorticity  streamlines  turning  toward  the  western  bound- 
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Figure  4.1:  Schematic  of  the  J-vector  fluxes  for  a  linear,  Stommel  and  Arons  solution  with 
Munk  boundary  layer. 


Figure  4.2:  J-vectors  for  run  Ul,  Re  =  5.  The  aspect  ratio  emphasizes  north-south  fluxes. 


Figure  4.3:  Streamfunction  of  J  for  run  Ul,  Re  =  5. 

ary  indicates  the  transfer  of  vorticity  to  that  boundary,  and  the  weak  poleward  and 
eastward  flux  in  the  interior  is  now  considerably  more  visible.  The  demarcation  line, 
xq  separating  J  having  these  two  destinies  is  not  drawn  in  the  figure,  but  its  existence 
and  location  is  easily  imagined. 


Weakly  Nonlinear 

Moving  to  the  next  stage  in  complexity  allows  nonlinearity  to  enter  the  dynamics, 
but  only  in  the  meridional  momentum  equation,  though  I  still  restrict  the  flow  to  be 
steady.  In  this  scenario  and  neglecting  the  meridional  flux  due  to  friction. 


Jed  = 
Jmn  ~ 

J  fric  = 


0 


-Ah 


(4.14a) 

(4.14b) 

(4.14c) 
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Now  the  total  advective  meridional  vorticity  flux  across  a  line  of  latitude  is 


iNarth  =  /  ^  +  1^) 

(4.15a) 

f^E  1  dv^ 

(4.15b) 

=  +  l  2  3^ -fa. 

=  yVK„,h  +  J  (xe)  -  (0)) 

(4.15c) 

If  the  boundary  conditions  are  no-slip,  this  reduces  to  the  same  vorticity  flux  as  in  the 
linear  case.  The  net  meridional  vorticity  flux  is  again  given  by  the  integrated  velocity 
advecting  planetary  vorticity  across  the  line  of  latitude. 

As  in  the  purely  linear  case,  the  equator  is  a  boundary  across  which  no  vorticity 
can  be  advected.  However,  one  feature  distinguishes  this  weakly  nonlinear  case  from 
the  purely  linear  one.  When  the  dynamics  are  linear,  there  is  zero  meridional  vorticity 
flux  at  all  longitudes  on  the  equator.  With  the  addition  of  weak  nonlinearity.  Equa¬ 
tion  4.14b,  there  can  be  cross-equatorial  flux  of  J  at  all  longitudes,  but  the  net  flux, 
when  integrated  across  the  domain,  must  be  zero. 

Consider  the  net  cross-equatorial  flux  of  vorticity  between  any  two  longitudes,  Xi 
and  X2- 

Isquator  =  V^{X2, 0)  -  0).  (4-16) 

Between  any  longitudes  of  identical  meridional  velocity,  there  is  no  net  flux.  Thus 
when  the  eastern  and  western  boundaries  are  used,  there  is  no  integrated  flux,  and 
similar  sub-domains  lie  between  all  locations  of  flow  reversal. 

The  streamfunction  of  J  for  a  weakly  nonlinear,  but  steady  numerical  experiment, 
run  U2  is  shown  in  Figure  4.4.  In  this  example,  the  streamlines  do  cross  the  equator 
near  the  intense  center  of  the  western  boundary  current,  but  return  both  to  the  west 
and  to  the  east  northward  such  that  the  integrated  flux  is  zero. 

Fully  Nonlinear 

Adding  further  complexity  by  allowing  nonlinearity  to  enter  the  dynamics  in  the  zonal 
momentum  equation  and  through  time  dependence  renders  the  problem  analytically 
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Figure  4.4:  Streamfunction  of  J  for  run  U2,  Re  =  20. 


intractable.  However,  it  is  possible  to  rely  on  the  numerical  calculations  in  this  limit 
Using  the  full  formulation  for  the  J-vectors  as  in  Equation  4.4a,  the  stream  function 
for  the  high  Reynolds  number  run  U8  is  shown  in  Figure  4.5  for  the  sub-domain  of 
the  system  bounded  by  latitudes  y  =  ±5  and  longitudes,  a:  =  0  and  x  =  2. 

The  immediately  striking  result  is  that  the  streamfunction  diagram  appears  re¬ 
markably  similar  to  that  in  the  weakly  nonlinear  run.  Although  the  system  has 
changed  quite  dramatically  from  the  linear  limit  (compare  again  Figures  3.7  and  3.10), 
the  dynamical  picture  given  simply  by  the  net  vorticity  flux  streamfunction,  is  not 
dramatically  different. 

Vorticity  flux  streamlines  enter  the  sub-domain  through  the  northern  boundary. 
Some  veer  westward  to  the  western  boundary,  and  others  diverge  to  the  east  and 
northward  around  the  mean  eddies  that  reside  in  these  regions.  Vorticity  transported 
along  these  paths  eventually  work  its  way  around  to  the  weak  interior  flow  and  out 
of  the  domain  across  latitude  y  —  5  (not  shown).  Some  streamlines  again  propagate 
across  the  equator,  with  one  extenting  quite  far.  However,  the  vorticity  flux  associated 
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Figure  4.5:  Streamfunction  of  J  for  run  U8,  Re  =  500. 


with  this  one  stream  tube  is  quite  small,  and  even  this  returns  northward  to  intersect 
the  western  boundary  very  near  to  the  equator.  Thus  even  in  the  highly  nonlinear 
case  examined,  vorticity  fluxes  southward  across  the  equator,  but  the  integrated  flux 
is  quite  small.  A  calculation  of  this  integrated  flux  gives  lEquator  =  0.002,  less  than  a 
percent  of  the  meridional  vorticity  flux  found  at  latitude  y  =  5. 

Although  the  linear  and  nonlinear  runs  appear  very  similar  using  this  diagnostic,  an 
important  distinction  exists  and  can  be  illustrated  by  decomposing  J  into  its  individual 
components  as  in  Equation  (4.5).  Figure  4.6  shows  the  contributions  to  J  by  each  of 
the  three  different  components  for  this  run  as  vector  plots.  Note  that  the  vector  plots 
are  not  to  the  same  scale,  but  can  be  used  to  infer  locations  where  different  processes 
are  active.  The  magnitudes  of  the  different  components  are  given  on  the  axes. 

In  Figure  4.6a,  the  mean  advection  of  mean  absolute  vorticity  shows  a  strong 
infliix  near  the  poleward  boundary  of  the  sub-domain.  Figure  4.6c  shows  the  flux 
due  to  friction.  Over  most  of  the  domain,  the  mean  frictional  flux  is  zero.  However, 
there  is  a  region  immediately  adjacent  to  the  western  boundary  where  this  flux  is 
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c)  Frictional 
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Figure  4.6:  J-vector  components,  and  Jjric,  for  run  U8,  Re  -  500.  Note  that  the 

quivers  are  not  to  same  scale,  but  the  maximum  of  each  vector  component  are  given  on  the  axes. 
The  extreme  aspect  ratio  causes  zonal  arrows  to  have  quite  small  heads  (particularly  noticeable  in 
the  frictional  plot). 
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significant.  It  is  difficult  to  see  this  in  the  figure  due  to  the  large  aspect  ratio  of  the 
plot  which  shrinks  the  heads  of  zonal  arrows  and  its  proximity  to  the  western  boundary. 
Careful  examination  reveals  that  particularly  in  the  northern  hemisphere,  though  also 
in  the  southern  hemisphere,  there  is  a  strong  westward  flux  of  vorticity  through  the 
western  boundary.  These  two  components  are  qualitatively  similar  to  those  of  the 
linear  and  weakly  nonlinear  dynamics  (not  shown).  However,  Figure  4.6b  shows  the 
contribution  by  the  eddy  field,  obviously  not  present  in  the  previous  experiments  since 
those  integrations  reached  steady-state  solutions. 

The  eddy  field  produces  a  diverging  vorticity  flux  north  of  the  equator.  Some 
vorticity  is  transferred  eastward,  where  it  merges  again  with  the  interior  flow  there. 
But  a  stronger  component  contributes  to  the  equatorward  and  westward  vorticity 
transport,  fluxing  vorticity  from  the  central  portion  of  the  boundary  current  southward 
and  westward,  where  it  can  be  removed  by  friction. 

Since  the  quiver  plots  are  not  to  the  same  scale  and  the  aspect  ratio  of  the  axes 
is  extreme,  a  better  reflection  of  the  relative  strengths  of  the  different  components  is 
given  by  the  divergences  of  the  individual  components  of  J.  Shown  in  Figure  4.7  is  the 
divergence  of  the  three  vector  fields  shown  in  Figure  4.6  in  a  still  narrower  sub-domain 
near  to  the  western  boundary.  A  positive  divergence  is  lightly  shaded  and  negative 
values  are  darkened.  The  divergence  of  the  different  components  shows  the  intensity 
of  J-vector  components  in  the  vicinity  of  the  equator  and  western  boundary  region. 
At  the  northern  portion  of  the  sub-domain,  the  convergence  of  meaji  J-flux  is  balanced 
by  a  frictional  divergence,  whereas  near  the  equator,  mean  and  frictional  fields  balance 
a  large  eddy-flux  convergence.  Slightly  further  from  the  boundary  region  is  a  strong 
eddy  divergence,  balanced  only  by  the  effect  of  the  mean  field. 

Although  one  effect  of  the  eddy-field  is  illustrated  by  the  J-vector  diagram  of  the 
previous  section,  additional  calculations  further  demonstrate  its  role  by  focussing  on 
the  contribution  of  the  eddy-field  to  the  net  westward  flux.  Since  excess  vorticity 
entering  the  sub-domain  through  the  northern  zonal  boundary  must  be  exported  out 
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a)  Mean 


c)  Frictional 
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X,  min(divJ)=  -1.3 


X,  min(divJ)=  -2.3 


X.  min(divJ)=-0.067 


Figure  4.7:  The  divergences  of  J-vector  components.  and  for  run  U8,  Re  =  500. 

Shading  is  identical  in  all  three  plots  and  maxima  are  given  on  the  axes.  Pure  white  corresponds  to 
values  less  than  -1.5. 


to  the  west,  it  is  interesting  to  inquire  which  processes  are  important  at  which  stages 
of  this  migration. 

Integral  measurements  are  determined  by  decomposing  Equation  4.6  into  its  four 
components: 


Jryn 

'  *^xue3t^y  ”1”  /  Jnorth^^ 

Va 

rvn  tXc 

I  Jcastdy  /  Jsouth^^ 

Jys  Jxw 

I  west  “h  Jnorth  “1“  least  I  south- 


(4.16b) 

(4.16c) 


where  small  letters  denote  boundaries  of  a  variable  subdomain  as  opposed  to  the 
full  basin  boundaries.  Integral  quantities  are  calculated  along  the  boundaries  Xyj  =  e, 
Xe  =  3,  y,  =  —5,  y„  =  5,  where  e  is  a  small  distance  from  the  western  boundary^.  The 
different  components  for  run  U8  gives  (rounded  to  hundredths) 


Inorth  —  —0.40 
Iwest  —  0.43 
Isouth  =  -0.03 

hast  =  0.00  . 


Noting  that  these  values  sum  to  zero  gives  some  measure  of  the  consistency  of  the 
model  in  accurately  carrying  vorticity.  As  expected,  the  flux  through  the  eastern  and 
southern  boundaries  of  the  subdomain  is  small  compared  to  the  northern  ajid  western 
ones.  If  the  subdomain  is  extended  to  the  eastern  boundary  there  is  a  slight  increa.se 
in  the  frictional  vorticity  flux  there  but  it  is  0(10“^),  consistent  with  expectations  of 
no  eastern  boundary  current. 

It  is  also  possible  to  decompose  each  geographic  I  component  into  its  mean,  eddy 
and  frictional  elements  and  observe  the  integrated  vorticity  budget  in  a  sub-domain. 
Table  4.1  displays  the  different  components  for  this  run.  As  expected  the  flux  through 

^  As  discussed  in  Appendix  A.4,  J-vectors  are  calculated  on  fc-points  of  the  model  and  are  therefore 
staggered  resting  just  off  the  boundary 
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3  tot 

Jmean 

I  eddy 

J/ric 

Iw 

.43 

0 

.01 

.42 

In 

-.40 

-.40 

0 

0 

Ie 

0 

-.01 

.01 

0 

Is 

O 

CO 

-.03 

0 

0 

Table  4.1:  Vorticity  budget  for  the  box  [0  3  -5  5]  for  run  U8. 

the  north  is  almost  entirely  due  to  the  mean  advection  of  mean  absolute  vorticity. 
And  the  flux  at  the  western  boundary  is  due  almost  entirely  to  friction. 

Although  there  appears  a  small  eddy  flux  at  the  western  boundary,  this  reflects 
the  small  distance,  a:  =  e,  at  which  the  measurement  is  made.  At  x  =  0  the  entire 
vorticity  flux  must  result  entirely  from  friction  since  the  mean  and  eddy  flux  at  a  solid 
boundary  is  exactly  zero. 

Though  useful,  this  calculation  only  gives  information  for  a  single  domain  and  is 
not  globally  illuminating.  By  scanning  the  western  edge  of  this  sub-domain  from  the 
western  boundary  eastward,  it  is  possible  to  explore  the  different  components  as  a 
function  of  longitude. 

A  scan  of  the  integrated  westward  flux,  for  northern  and  southern  hemi¬ 

spheres  of  this  run  is  shown  in  Figure  4.8.  In  solid  is  the  total  westward  flux  of 
vorticity,  dashed  is  that  due  to  friction.  The  effect  of  the  mean  field  is  the  dotted 
curve  and  the  dot-dcished  is  the  westward  flux  due  to  the  eddies.  The  heavy  solid 
line  shows  the  residual  of  the  sum  of  the  different  integrated  components  around  the 
whole  sub-domain,  F.  In  both  figures,  this  residual  is  quite  close  to  zero  for  all  longi¬ 
tudes,  consistent  with  Equation  4.6.  In  both  hemsipheres,  there  exists  a  very  intense 
westward  flux  due  to  friction  at  the  western  boundary.  However,  the  frictional  bound- 
ary  layer  is  narrow  compared  to  the  width  of  the  boundary  current  itself  (compare 
Slf  =  0.044  to  6 1  =  0.1)  in  this  high  Reynolds  number  run.  With  a  large  portion  of  the 
full  boundary  layer  passing  outside  the  region  where  friction  can  act,  the  system  ad¬ 
justs  by  creating  an  eddy  field,  to  the  east  of  the  frictional  boundary  layer,  to  transport 
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vorticity  from  the  interior  of  the  boundary  current  to  the  frictional  sublayer.  This  role 
is  illustrated  by  the  strong  hump  in  the  dot-dashed  curve  at  x  =  0.1.  Interestingly, 
in  the  northern  hemisphere  where  the  vorticity  constraint  in  the  linear  and  weakly 
nonlinear  systems  is  most  compelling,  the  mean  field  acts  in  fact  to  advect  vorticity 
away  from  the  boundary.  The  eddy  and  frictional  components  must  compensate  not 
only  for  the  narrowing  of  the  frictional  sublayer,  but  also  for  this  competing  effect  of 
the  mean  field.  In  the  southern  hemisphere,  the  mean  field  flux  has  a  different  sign 
due  to  the  change  in  the  sign  of  the  Coriolis  parameter,  but  the  effect  of  the  eddy 
field  is  still  quite  apparent  in  the  narrow  region  immediately  adjacent  to  the  frictional 
sublayer. 

To  examine  the  role  of  the  eddies  in  transporting  vorticity  meridionally,  the  inte¬ 
grated  flux  of  vorticity  between  x  =  0  and  a:  =  3  as  a  function  of  latitude  is  shown 
in  Figure  4.9.  There  is  a  large  meridional  vorticity  flux  at  high  latitudes  due  almost 
entirely  to  the  mean  field,  but  this  flux  decays  steadily  with  latitude  and  is  virtually 
zero  at  the  equator.  A  small  residual  remains  due  to  the  eddy  flux.  This  figure  demon¬ 
strates  that  in  the  mean,  zonally  integrated  sense,  a  system  constrained  to  transport 
fluid  across  the  equator  successfully  changes  the  sign  of  the  vorticity  of  the  fluid  in 
the  boundary  current,  even  with  a  very  small  magnitude  of  dissipation.  By  the  time 
the  fluid  reaches  the  equator,  the  large  anomalous  absolute  vorticity  of  the  boundary 
current  is  largely  reduced  to  a  small,  equatorial  value. 

The  broader  effect  of  varying  the  Reynolds  number  can  be  determined  by  examining 
the  ratio  of  the  peak  of  the  dashed  curve  to  the  maximum  of  the  integrated  frictional 
flux  in  Figure  4.8.  One  difficulty  in  defining  this  index  is  that  it  depends  on  the  choice 
of  the  maximum  latitude,  t/„,  for  integration.  The  frictional  flux  extends  all  along 
the  western  boundary,  whereas  the  maximum  westward  eddy  flux  is  greatest  within  a 
few  deformation  radii  of  the  equator  (see,  for  example,  the  lightly  shaded  regions  in 
Figure  4.7).  Interestingly,  the  extent  of  the  region  of  enhanced  eddy  flux  is  similar  for 
a  range  of  Reynolds  numbers,  and  lies  between  the  equator  and  a  latitude  between 
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Figure  4.9:  — for  run  U8,  Re  =  500.  Lines  are  as  in  Figure  4.8.  The  large  meridional  flux  of 

vorticity  at  latitude  y  =  5  decreases  to  nearly  zero  at  y  =  0.  Note  that  the  integral  in  this  calculation 
extends  across  the  basin. 


y  =  3  and  y  =  4  for  this  model  configuration. 

For  each  Reynolds  number,  I  estimate  the  maximum  latitude  for  integration  as 
that  location  where  the  integral  of  the  westward  eddy  flux  has  reached  98%  of  its 
value  at  y  =  5,  which  is  assumed  to  be  close  to  its  maximum  value,  and  define  the 
index, 


(4.14) 


to  represent  an  integrated  and  normalized  measure  of  the  strength  of  the  eddy  field 
in  transporting  vorticity  to  the  boundary.  This  index  is  calculated  from  numerical 
experiments  and  plotted  in  Figure  4.10  as  a  function  of  Reynolds  numbers.  At  low 
Re,  the  integrated  eddy  flux  is  quite  small.  The  star  corresponding  to  Re  =  32  (run 
U2.A)  has  an  index,  E  =  0.003,  and  indicates  that  the  critical  Reynolds  number  for 
instability  to  occur  is  slightly  less  than  this  value.  As  the  Reynolds  number  increases, 
the  integrated  effect  of  the  eddies  increases  also,  showing  the  growing  role  of  the  eddies 
in  transporting  vorticity  westward.  Included  in  this  plot  are  two  runs  at  Re  =  150 
which  examined  variation  in  basin  dimension.  Results  from  both  run  U5.X  (with 
(Lx,  Ay  =  (15,20))  and  run  U5.Y  (with  (Lx,  Ay)  =  (10,30))  show  remarkably  good 
agreement  with  one  another  and  in  the  trend  of  this  integral  calculation. 
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Figure  4.10:  Index  E,  defined  by  Equation  4.14,  showing  the  ratio  of  maximum  integrated 

westward  eddy  flux  of  vorticity  to  the  integrated  dissipative  flux  as  a  function  of  Reynolds  number. 
Integration  is  extended  from  the  equator  to  the  latitude  where  the  westward  eddy  flux  reaches  98%  of 
its  maximum  value.  The  circle  corresponds  to  a  run  with  non-uniform  upwelling  (UFS).  Runs  U2.A 
U4,  U5.X,  U5.Y,  U6,  U7,  U8,  U9,  and  N2  are  included. 

The  figure  suggests  that  the  integrated  eddy  flux  increases  approximately  expo¬ 
nentially  from  zero  at  a  critical  value  of  Rcc  ~  30  to  a  value  of  more  than  1  at 
Re  =  500.  This  critical  value  is  determined  empirically,  but  is  not  drastically  different 
from  other  calculations  using,  for  example,  a  nonuniform  upwelling  distribution.  The 
approximate  exponential  growth  of  the  normalized  flux  is  not  predicted  by  a  theory. 

4.1.3  Balances  from  the  Time-Mean  Vorticity  Equation 

The  vorticity  balance  of  the  linear  system  is  obvious  from  the  conditions  of  linearity. 
Specifically,  planetary  vorticity  advection  is  balanced  by  vorticity  diffusion.  As  the 
system  becomes  more  nonlinear,  the  boundary  layer  structure  is  not  obvious,  but  can 
be  diagnosed  from  the  model.  This  section  examines  the  structure  of  the  boundary 
layer  as  a  function  of  latitude  in  the  high  Reynolds  number  runs. 
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The  time-averaged  vorticity  equation  (4.1)  is  written  again, 


u  •  VC  +  v  +  (C  +  y)  V  •  u  =  AhVK  (4.15) 

where  the  different  terms  correspond  to  relative  vorticity  advection  (RVA),  planetary 
vorticity  advection  (PVA),  stretching  (STR),  and  friction  (FRIG).  The  eddy  terms 
have  not  been  separated  from  the  different  components  of  vorticity  as  the  purpose  of 
this  section  is  to  determine  the  time-mean  structure  of  the  western  boundary  current. 
The  numerical  diagnosis  of  these  terms  according  to  the  discretization  of  the  model  is 
discussed  in  Appendix  A.3. 

For  very  large  Reynolds  number,  the  inertial  boundary  layer  scale,  8i  =  y/U ,  is 
larger  than  the  Munk  scale.  In  this  case  the  system  divides  into  multiple  regions  with 
different  balances  presiding.  The  different  components  of  the  time-averaged  vortic¬ 
ity  equation  for  a  subdomain  of  run  U7  near  the  western  boundary  are  shaded  and 
contoured  in  Figure  4.11.  Frictional  dissipation  (Figure  4. lid)  occurs  entirely  in  the 
intense  frictional  boundary  layer  that  resides  immediately  adjacent  to  the  western 
boundary.  Its  magnitude  and,  to  some  extent,  zonal  scale  decrease  with  latitude. 
The  planetary  vorticity  advection  ((Figure  4.11b)  is  proportional  to  v  and  therefore 
reveals  the  zonal  width  of  the  full  boundary  layer.  The  boundary  current  in  this  non¬ 
linear  run  is  relatively  narrow  in  the  northern  portion  of  the  sub-domain,  but  widens 
considerably  near  the  equator.  In  addition,  the  magnitude  of  the  planetary  vorticity 
advection  decreases  substantially  by  this  latitude.  The  relative  vorticity  advection 
((Figure  4.11a)  is  characterized  by  an  intense  layer  adjacent  to  the  western  boundary 
but  with  an  oppositely  signed,  weaker  contribution  to  its  east.  The  effect  of  stretching 
is  found  in  Figure  4.11c  and  is  surprisingly  small  in  all  regions  of  this  sub-domain. 
Of  course,  the  contribution  of  stretching  to  the  vorticity  balance  is  significant  in  the 
interior  (roughly  2  to  3  deformation  radii  from  x  =  0)  and  away  from  the  equator, 
where  stretching  balances  planetary  vorticity  advection. 

To  better  understand  the  balances  from  this  plan  view,  the  zonal  profiles  of  the 
different  components  at  4  latitudes  are  shown  in  Figure  4.12.  The  upper  panel  shows 
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the  different  components  at  y  =  5  and  reveals  a  split  of  the  boundary  current  into  two 
dynamical  regimes.  The  inner  layer  shows  a  three-way  balance  of  friction  (dotted), 
relative  vorticity  advection  (solid),  and  planetary  vorticity  advection  (dashed).  There 
is  also  an  outer  layer  that  is  primarily  an  inertial  balance  of  RVA  and  PVA.  As  the 
latitude  decreases,  the  inner  layer  narrows  considerably,  the  outer  layer  widens,  and 
the  effect  of  planetary  vorticity  advection  in  the  inner  layer  decreases.  This  balance 
remains  valid  far  south  of  the  equator  as  well. 

From  the  input  parameters  it  is  possible  to  estimate  the  expected  scales  of  the 
boundary  layers.  For  run  U7,  the  Munk  scale  is  6m  =  0.0667,  the  inertial  scale, 
assuming  a  Stommel  and  Arons  interior  flow  as  in  Equation  2.19,  is  0.1.  However, 
the  three-way  balance  of  Figure  4.12a  suggests  a  scaling  of  the  system  such  that 
=  6m-  Thus  in  this  northern  portion  of  the  subdomain,  there  exists  a  slightly 
narrower  inertial  current  than  expected  from  Equation  2.19.  Moving  further  south  in 
the  basin,  the  two-way  balance  of  RVA  and  FRIG  suggests  an  inertio-frictional  balance 
in  the  inner  layer: 


Using  6lf  =  0.03  taken  approximately  from  the  figure.  Equation  4.16  suggests  6i  ~  0.3. 
This  value  compares  sensibly  with  the  scale  of  the  outer  boundary  layer  in  Fig¬ 
ures  4.12b,  c,  and  d,  using  the  zero  crossing  of  the  planetary  advection  (double-dashed) 
curve  at  s  =  0.6  as  a  reference  location.  Again  this  balance  differs  from  that  suggested 
by  linear  theory. 

The  vorticity  balance  suggests  a  widening  of  the  inertial  layer  with  decreasing 
latitude.  As  the  parcels  approach  the  equator,  they  are  swept  out  into  the  interior 
by  the  eddy  that  resides  there  in  the  mean.  The  effect  of  this  eddy  is  observed  in 
the  widening  of  the  scale  of  the  PVA  term  in  Figure  4.12b,  c,  and  d,  and  by  the  very 
rough  scaling  presented  above. 
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Figure  4.12:  Zonal  profiles  of  the  different  components  of  the  vorticity  equation,  shown  in  Fig¬ 
ure  4.11  at  4  different  latitudes,  (a)  y  =  5,  (b)  2.5,  (c)  0,  and  (d)  -2.5. 
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4.2  Localized  Up  welling 


A  critical  aspect  of  cross-equatorial  flow  is  the  determination  of  the  actual  cross- 
equatorial  mass-fliix.  In  this  work,  it  is  specified  through  the  upwelling  distribution, 
lu*.  The  results  presented  thus  far  concern  a  very  specific  system  in  which  the  upwelling 
distribution  is  constant,  and  therefore  may  be  considered  a  potential  weakness  of  this 
study.  Indeed,  although  Stommel  and  Arons  (1960a)  considered  a  general  form  for  the 
upwelling  and  directly  examined  both  uniform  and  non-uniform  distributions  in  their 
paper,  the  work  is  frequently  referred  to  as  pertaining  to  uniform  upwelling  only. 

In  an  effort  to  demonstrate  that  the  results  presented  here  are  robust  to  different 
upwelling  distributions,  I  conduct  an  additional  experiment  that  considers  the  most 
diametrically  opposed  forcing,  a  localized  tn*  (Equation  2.15).  In  this  configuration, 
all  the  fluid  entering  in  the  northern  hemisphere  is  destined  to  exit  the  basin  through 
a  sink  entirely  in  the  southern  hemisphere.  Thus  all  of  the  fluid  entering  the  system 
is  constrained  to  cross  the  equator.  One  limitation  of  this  configuration  is  that  the 
volume  flux  of  fluid  out  of  the  sink  must  not  be  so  large  that  the  layer  depth  vanishes. 
Although  this  is  not  a  problem  at  low  Reynolds  number,  it  can  be  significant  as  the 
volume  of  fluid  transported  through  the  system  increases.  As  a  result,  the  sink  has  a 
larger  cross-sectional  area  than  the  source  (10  x  10  versus  5  X  10  deformation  radii), 
and  iw*  is  positive  for  all  latitudes  south  of  y  =  -5.  A  plot  showing  the  upwelling 
distribution  is  presented  in  Figure  4.13. 

The  system  evolves  much  as  the  system  with  distributed  upwelling.  There  is  an 
initial  creation  of  a  mound  of  fluid  in  the  forcing  region  with  subsequent  Kelvin  and 
Rossby  waves  setting  up  the  interior  circulation.  Unlike  the  previous  runs,  the  inte¬ 
rior  circulation  is  now  absent,  but  there  remains  after  the  spin- up  a  residual  western 
boundary  current  that  carries  fluid  from  source  region  to  sink. 

As  in  the  uniform  upwelling  ca.se,  the  system  becomes  unstable  at  higher  Reynolds 
number,  and  eddies  develop  in  the  vicinity  of  the  equator.  In  Figure  4.14  is  shown  the 
instantaneous  velocity  and  height  calculations  from  a  Re  =  75  run.  The  accompanying 
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Figure  4.13:  The  upwelling  distribution  in  the  basin  for  run  N2.  The  source  is  located  in  the 
northern  hemisphere  and  the  sink  in  the  southern,  with  zero  vertical  mass  flux  in  between. 

potential  vorticity  field  is  presented  in  Figure  4.15.  Eddies  form  immediately  north  of 
the  equator,  propagate  steadily  to  the  south,  and  dissipate  vorticity  along  their  path. 
In  the  mean,  the  boundary  current  widens  at  the  equator  (Figure  4.16).  As  before, 
this  feature  is  associated  with  a  recirculation,  though  displaced  somewhat  south  of  its 
location  in  the  uniform  upwelling  case.  Again,  this  recirculation  is  not  indicative  of 
trapped  fluid,  but  rather  of  a  location  of  enhanced  eddy  activity. 

The  vorticity  flux  streamfunction,  defined  by  Equation  4.13  takes  a  similar 
form  as  in  the  uniform  upwelling  case  (Figure  4.17).  Consistent  with  the  requirement 
that  all  of  the  fluid  passing  through  the  system  must  cross  the  equator  and  therefore 
change  its  vorticity,  most  of  the  vorticity  entering  the  subdomain  of  the  basin  through 
the  northern  boundary  exits  through  the  west  (a  small  portion  is  fluxed  back  to  the 
north  in  the  counter-current  that  resides  just  east  of  the  western  boundary  current). 
Calculating  the  index  E  for  this  run  (integrating  only  to  y  =  1.75  where  the  integrated 
eddy  flux  is  nearly  constant)  gives  a  value  of  E  =  0.33,  which  is  slightly  lower,  but 
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Figure  4.17:  The  vorticity  flux  streamfunction  for  run  N2.  Unlike  the  uniform  upwelling  runs, 
most  of  the  vorticity  entering  the  subdomain  of  the  basin  at  y  =  5  leaves  through  the  western 
boundary. 

reasonably  consistent  with  the  behavior  in  Figure  4.10  and.  shown  as  the  circle  in  that 
plot. 

The  details  of  the  localized  upwelling  run  are  somewhat  different  from  those  of 
the  previous  runs.  However,  the  basic  mechanism  is  identical.  The  ratio  of  inertial 
to  frictional  boundary  current  increases  with  the  Reynolds  number.  At  a  critical 
number,  approximately  Re  =  30,  the  system  goes  unstable,  and  the  eddies  provide 
the  mechanism  for  the  necessary  lateral  flux  of  vorticity  to  enable  the  conversion 
necessary  for  cross-equatorial  flow.  The  upwelling  specifies  the  mass  that  must  flux 
across  the  equator.  This  section  shows  that  regardless  of  the  distribution  of  upwelling, 
there  exists  a  sufficient  eddy  flux  to  allow  that  transport. 
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4.3  Pree-slip  Boundary  Conditions 


Although  float  data  in  the  equatorial  Atlantic  ocean  indicate  that  the  velocity  par¬ 
allel  to  the  coast  becomes  small  at  the  coast  itself  (Figure  1.3),  I  now  consider  the 
sensitivity  of  the  results  to  this  quantity  by  analyzing  how  the  circulation  changes 
with  a  different  boundary  condition.  Specifically,  I  examine  the  the  effect  of  a  free- 
slip  boundary.  As  can  be  expected  from  other  studies  (e.g.,  Marshall  (1981)  or  the 
discussion  in  lerley  (1990)),  this  change  alters  the  general  circulation,  the  boundary 
layer  structure,  and  the  eddy  fluxes  from  those  using  the  no-slip  condition.  This  sec¬ 
tion  demonstrates  that  although  some  features  do  change  with  boundary  conditions, 
the  physical  mechanism  of  eddy  flux  transport  remains  essential  to  the  dynamics  of 
nonlinear  cross-equatorial  flow. 

Before  considering  the  fully  nonlinear  case,  it  is  instructive  to  revisit  the  steady, 
weakly  nonlinear  discussion  of  Section  4.1.2.  The  fundamental  dynamical  change  to  the 
system  due  to  the  new  boundary  condition  is  apparent  in  the  integrated  meridional 
vorticity  flux  (Equation  4.15c).  With  free-slip  boundary  conditions,  the  constraint 
that  u^(x)  vanish  at  the  western  boundary  no  longer  applies,  and  the  net  meridional 
vorticity  flux  can  be  conSderable,  even  at  the  equator.  In  turn,  the  relaxation  of  this 
constraint  allows  the  system  a  greater  meridional  span  over  which  to  remove  vorticity 
from  the  basin,  and  therefore  the  local  westward  vorticity  flux  is  reduced  from  the 
no-slip  system. 

Consider  the  Re  =  200  run  UFS  whose  upwelling  is  uniform,  but  the  free  slip 
condition  is  applied.  Mean  velocity  and  height  fields  are  shown  in  Figure  4.18.  The 
mid-latitude  interior  presents  the  familiar  Stommel  and  Arons  (1960a)  eaistward  and 
poleward  flow.  Near  the  western  boundary  rests  an  intense  southward  flow  on  and 
immediately  adjacent  to  the  boundary,  with  an  intense  counterflow  just  to  its  east. 
The  mean  recirculation  that  had  been  present  in  the  vicinity  of  the  equator  in  the 
no-slip  runs  now  rests  quite  far  south  in  the  southern  hemisphere  and  is  considerably 
elongated  meridionally,  compressed  zonally. 
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Figure  4.19:  The  vorticity  flux  streamfunction  for  a  sub-domain  of  run  UFS. 

The  vorticity  flux  streamfunction,  shows  a  large  southward  vorticity  flux  at  the 
northern  end  of  the  subdomain,  as  before  (Figure  4.19).  However,  now  the  stream¬ 
lines  pass  smoothly  across  y  =  0,  eventually  turning  westward  out  of  the  basin,  but 
extending  south  beyond  the  equator  several  deformation  radii. 

Two  integrated  flux  measurements  are  presented  in  Figure  4.20.  In  the  upper  panel 
is  the  net  meridional  vorticity  flux,  considerable  at  the  northern  edge  of  the  domain 
and  decreasing  with  latitude.  This  quantity  vanishes  at  y  =  —6.  In  the  lower  panel 
is  the  net  westward  flux  between  y  =  —6  and  y  =  6.  Near  to  the  boundary  the  total 
westward  flux  is  primarily  due  to  friction.  To  be  sure,  there  is  a  small  contribution 
to  the  westward  flux  in  a  turbulent  boundary  layer,  visible  by  the  small  hump  in 
the  dot-dashed  curve  near  x  =  0.1,  This  flux  is  strongest  between  latitudes  y  =  — 4 
and  y  =  —  7,  far  south  of  the  equator.  Relaxing  the  constraint  of  no  cross-equatorial 
meridional  flux  allows  friction  to  dominate  the  boundary  layer  physics  except  in  the 
vicinity  of  y  =  — 6  where  the  zonally  integrated  vorticity  flux  is  zero. 

The  significant  role  that  friction  plays  in  the  boundary  layer  dynamics  is  also  found 
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Figure  4.20:  The  integrated  westward  and  northward  vorticity  fluxes  for  run  UFS.  The  dashed 
curve  corresponds  to  the  frictional  flux,  the  dotted  to  mean  advection,  the  dot-dashed  to  the  eddy 
flux,  and  the  light  solid  curve  is  their  sum.  a)  The  integral  extends  from  y  =  — 6  to  y  =  5,  and 
the  heavy  solid  curve  is  the  sum  of  net  fluxes  through  a  sub-domain  bounded  by  these  latitudes, 
the  longitude  in  the  figure,  and  x  =  2.  b)  The  integral  extends  zonally  across  the  bzisin,  and  the 
sub-domain  extends  meridionally  from  the  latitude  in  the  figure  to  y  =  5. 
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in  the  boundary  layer  structure.  Shown  in  Figure  4.21  are  the  different  components 
of  the  vorticity  equation  near  the  western  boundary  at  4  latitudes.  In  the  northern 
portion  of  the  basin  at  7/  =  5,  the  boundary  layer  consists  of  two  layers,  an  inner 
viscous  sublayer  of  the  inertio-frictional  scale  {8lj  =  0.06)  and  an  outer  inertial  layer. 
The  meridional  flow  decays  eastward  from  the  western  boundary  and  changes  sign 
near  x  =  0.35.  Near  the  equator,  the  full  boundary  layer  width  has  decreased  and  the 
meridional  velocity  changes  sign  near  x  =  0.2,  and  the  inner  inertial  layer  has  been 
subsumed  into  a  wider  frictional  layer.  This  frictional  balance  extends  south  of  the 
equator,  and  is  shown  for  example  at  y  =  —2.5.  Far  south  of  the  equator  at  y  =  —6, 
the  boundary  layer  remains  narrow,  but  the  balance  is  more  of  a  three  term  balance, 
both  in  the  inner  and  outer  portions  of  the  southward  flowing  current.  It  is  at  this 
latitude  that  the  structure  of  the  boundary  current  most  resembles  the  near  equatorial 
balances  of  the  no-slip  case. 

Although  the  transfer  of  vorticity  to  the  boundary  is  largely  achieved  by  the  fric¬ 
tional  boundary  layer,  eddies  do  become  significant  far  south  of  the  equator,  and  also 
in  the  interior.  This  interior  process  is  represented  in  the  strong  peaks  of  integrated 
eddy  flux  in  Figure  4.20b,  eastward  at  x  =  0.35  arid  westward  at  x  =  0.9.  Both 
functions  of  the  eddy  field  are  further  highlighted  in  Figure  4.22  which  shows  the  com¬ 
ponents  of  the  J- vectors  due  to  the  eddy  field.  Immediately  adjacent  to  the  boundary, 
the  eddy  field  flux  is  negligible,  except  in  a  small  region  near  the  southern  portion  of 
the  domain.  Furthermore,  the  interior  flux  is  easily  visible  in  the  strong  convergence 
of  vorticity  near  x  =  0.6.  This  convergence  demarcates  the  eddy  flux  transfering  vor¬ 
ticity  to  the  northward-flowing  counter-current  from  the  western  boundary  current. 
Both  the  boundary  current  and  the  counter-current  cross  contours  of  mean  potential 
vorticity.  Associated  with  this  crossing  is  a  necessary  vorticity  flux  divergence,  pre¬ 
dominantly  due  to  the  eddy  field.  Thus  as  in  the  no-slip  case,  eddies  move  vorticity 
to  the  boundary  (in  a  small  region  in  the  southern  hemisphere),  and  they  are  also 
responsible  for  the  interior  transfer  from  meridional  flow  to  meridional  flow  and  drive 
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Figure  4.21:  Different  components  of  the  vorticity  equation  at  4  latitudes  in  run  UFS:  (a)  v  =  5 
(b)  y  =  0,  (c)  j/  =  —2.5,  and  (d)  y  =  —6. 


Figure  4.22:  Eddy  flux  vectors  for  run  UFS. 

the  recirculation. 

It  is  clear  from  observations  such  as  Figure  1.3  that  the  free  slip  boundary  condition 
does  not  reflect  nature;  however,  setting  u  =  0  exactly  on  the  boundary  may  also  be 
inadequate,  given  the  coarse  resolution  of  even  this  model.  It  is  important  to  draw 
from  the  above  example  the  fact  that  although  some  features  of  the  circulation  have 
changed,  fundamental  aspects  are  similar.  No  longer  constrained  to  flux  vorticity 
out  of  the  basin  prior  to  reaching  the  equator,  the  system  can  maintain  a  frictional 
boundary  layer  over  a  meridionally  extended  portion  of  the  basin,  without  invoking 
the  eddy  field.  However,  ultimately,  the  system  confronts  a  similar  constraint  as  in 
the  no-slip  case.  At  some  latitude,  the  advective  vorticity  flux  of  the  steady  system 
must  vanish.  This  special  latitude  is  given  by 

y  I  vdx  =  -^^(0)  (4-17) 

Jo  2 

which  must  be  internally  determined  by  the  system.  In  the  above  example,  this 
equality  is  satisfied  at  approximately  y  =  —5.1.  It  is  this  constraint  which  the  steady 
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system  (frictional  by  definition)  can  not  accomodate,  and  therefore  requires  the  eddy 
field  in  this  region  to  aid  in  the  vorticity  transfer.  The  boundary  condition  simply 
determines  that  latitude  where  the  net  meridional  flux  of  vorticity  must  approximately 
vanish.  For  no-slip,  it  resides  at  the  equator.  For  free-slip  and  strongly  forced,  it  can 
be  some  distance  south.  Naturally,  an  intermediate  boundary  condition  will  produce 
an  intermediate  critical  latitude,  and  the  growth  of  an  eddy  flux  at  that  location. 

4.4  Barotropic  Limit 

One  surprising  result  above  is  that  the  stretching  term  in  the  mean  vorticity  equa¬ 
tion  of  Section  4.1.3  appears  inconsequential  in  the  dynamical  balance  of  the  western 
boundary  current  (refer  again  to  Figure  4.12).  The  stretching  term  is  an  essential 
element  of  the  shallow-water  system,  the  critical  feature  that  distinguishes  it  from 
purely  barotropic  models  of  the  ocean.  It  is  also  the  stretching  that  distinguishes  the 
equator  in  the  shallow-water  model  from  all  other  latitudes.  This  fact  can  be  seen 
from  the  inviscid  quasi-geostrophic  vorticity  equation  in  which  the  stretching  term  is 
absent. 

Ct  +  u-VC  +  v  =  0.  (4.18) 

Without  /  in  the  governing  equation,  the  QG  system  considers  only  variations  in  /, 
not  its  absolute  value.  Returning  to  the  vorticity  calculations,  the  stretching  term 
appears  to  be  unimportant,  but  can  not  be  neglected  or  the  dynamics  would  behave 
quite  differently. 

A  demonstration  of  this  fact  is  found  by  examining  the  barotropic  limit  of  the 
model,  in  which  the  boundary  layer  is  small  compared  to  the  deformation  radius. 
Figure  4.23  shows  the  instantaneous  height  and  velocity  field  for  model  run  UBar 
i,R^  100)  at  t  =  500(/5Z/£))  Note  that  a  quasi-steady  state  is  achieved  more 
quickly  in  this  run  because  of  the  fast  barotropic  Rossby  waves  that  return  westward 
across  the  basin.  The  figure  is  dominated  by  long  basin  modes  propagating  westward 
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in  both  hemispheres.  However,  in  the  vicinity  of  the  western  boundary  current  are 
again  eddy  features.  However,  unlike  the  baroclinic  case,  the  variability  develops 
immediately  south  of  the  forcing  region,  and  extends  steadily  along  the  boundary  at 
least  beyond  y  ^  —0.2. 

Thus  when  the  boundary  current  is  small  compared  to  the  deformation  radius,  the 
system  is  basically  barotropic,  and  the  model  can  not  distinguish  any  unique  equatorial 
physics.  In  the  more  baroclinic  runs  described  earlier  in  the  chapter,  the  equator 
clearly  does  play  a  role.  Surprisingly,  from  the  time-mean  calculations,  thickness 
advection  appears  to  be  small  in  the  mean  vorticity  balance  of  the  western  boundary 
current.  Of  course,  stretching  is  important  to  the  interior  mid-latitude  balance,  but 
this  reflects  the  linear  effect  of  the  forcing,  and  not  the  maintenence  of  the  height  field 
fundamental  to  the  shallow-water  model. 

It  is  possible  that  the  the  system  recognizes  the  equator  through  the  time-dependent 
stretching  term,  which  vanishes  in  the  time  mean.  Indeed  it  is  found  in  Chapter 
5  that  this  term  is  essential  to  the  development  of  eddies  preferrentially  in  the  vicinity 
of  the  equator.  In  the  time-dependent  model,  this  form  of  stretching  could  contribute 
to  the  vorticity  balance  at  high  latitudes,  where  /  is  large,  and  the  system  is  most 
baroclinic.  Near  the  equator,  /  is  weak,  and  the  effect  of  relative  vorticity  advection 
may  be  more  dominant.  However,  cursory  examination  of  instantaneous  fields  in  the 
strongly  nonlinear.  Re  =  500  run,  shows  little  correlation  between  the  time  evolution 
of  vorticity  and  this  stretching  term. 

It  is  important  to  remember  that  the  Kelvin  waves,  so  vital  in  the  development  of 
the  circulation,  require  the  presence  of  stretching  term.  From  the  strong  variability 
that  exists  in  the  equatorial  zone  in  Figure  3.10,  it  appears  that  these  waves,  and 
therefore  the  stretching  term,  may  also  play  an  important  role  contributing  to  the 
quasi-steady  state  as  well,  though  I  have  not  analyzed  this  aspect  of  the  model.  It  re¬ 
mains  an  interesting  question  how  exactly  the  change  in  sign  of  the  Coriolis  parameter 
influences  the  dynamics  of  the  vorticity  balance. 
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Figure  4.23:  Instantaneous  height  contours  and  velocity  quivers  for  the  barotropic  model  run 
UBar  {Re  =  100)  at  f  =  500.  Notice  that  the  variability  in  the  western  boundary  current  forms  just 
south  of  the  intense  forcing  region  in  the  northwestern  corner. 


4.5  Summary 


Model  integrations  of  Chapter  3  provide  a  qualitative  description  of  the  behavior  of 
fluid  in  cross-equatorial  flow  under  a  range  of  conditions.  The  distribution  of  up  welling 
sets  the  volume  flux  of  fluid  that  must  cross  the  equator,  and  the  system  responds 
by  creating  an  eddy  field  as  the  ratio  of  the  volume  flux  of  fluid  to  the  viscosity 
increases.  The  formalism  of  the  J- vector  analysis  introduced  in  this  chapter  reveals 
and  quantifies  the  dynamical  role  of  the  eddy  field  that  appears  in  the  numerical 
experiments.  The  eddies  serve  to  transfer  vorticity  between  inertial  and  frictional 
portions  of  the  boundary  current,  enabling  the  potential  vorticity  transformation  that 
permits  the  flow.  The  strength  of  the  turbulent  boundary  layer  increases  with  the 
Reynolds  number  perhaps  without  bound  and  at  least  to  Re  =  1000.  The  behavior 
of  the  system  is  quite  robust  to  different  up  welling  configurations,  and  though  the 
circulation  changes  considerably  using  free-slip  boundary  conditions,  the  eddy  field 
again  plays  a  dynamical  role  in  its  vorticity  budget.  In  the  vicinity  of  the  equator, 
the  structure  of  the  no-slip  uniform  upwelling  boundary  current  is  divided  into  two 
portions,  an  outer  inertial  boundary  layer,  and  an  inner  layer  obeying  an  inertio- 
viscous  balance. 

The  analysis  above  reflects  its  Eulerian  perspective.  It  states  how  vorticity  moves 
in  the  frame  of  reference  of  the  b2tsin.  Speculating  on  the  Lagrangian  behavior  of 
individual  fluid  parcels  provides  an  alternative  description  of  cross-equatorial  flow. 
In  strongly  nonlinear  flow,  a  parcel  travels  southward  within  the  western  boundary 
current,  its  vorticity  is  partially  conserved,  only  weakly  modified  through  a  dissipative 
process.  Near  the  equator,  this  parcel  is  swept  into  an  intense  cyclonic  structure  and 
immersed  in  a  sea  of  eddies  having  both  positive  and  negative  senses  of  circulation. 
As  the  eddy  sweeps  past  these  eddies,  and  importantly,  also  past  the  boundary,  it 
inevitably  passes  through  narrow  zones  of  intense  velocity  shear,  where  the  vorticity 
is  transferred,  from  parcel  to  parcel,  and  ultimately,  from  parcel  to  boundary.  The  path 
of  any  given  parcel  is  extremely  complicated  and  difficult  to  model,  and  the  trajectories 
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of  neighboring  parcels  diverge  rapidly  once  they  enter  the  turbulent  region.  Once  the 
parcel  has  transferred  all  of  its  anomalous  vorticity,  it  can  merge  smoothly  with  a 
relatively  weak  interior  flow,  either  in  the  northern  or  southern  hemisphere,  where  the 
dynamics  are  more  linear. 
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Chapter  5 


The  Tropically  Enhanced 
Instability  of  the  Western 
Boundary  Current 

5.1  Introduction 

Central  to  the  discussion  of  cross-equatorial  flow  in  the  deep  ocean  is  the  development 
of  the  eddy  field  that  arises  as  the  system  becomes  increasingly  nonlinear.  Associated 
with  these  eddies  is  the  necessary  transfer  of  vorticity  between  different  parts  of  the 
boundary  current,  which  allows  fluid  parcels  to  modify  their  potential  vorticity  along 
their  meridional  journey.  The  existence  of  these  eddies  and  their  dynamical  role  are 
now  well-established  by  the  preceding  experiments,  but  their  cause  is  not  explained. 
The  goal  of  this  chapter  is  to  analyze  the  instability  of  the  tropical  western  boundary 
current  and  explain  the  mechanism  responsible  for  the  eddy  generation. 

Strong  eddy  activity  in  the  tropical  oceans  is  not  purely  a  numerical  artifact  but  is 
also  found  in  nature.  As  mentioned  in  the  Introduction  to  the  thesis,  many  observa¬ 
tions  indicate  that  the  equatorial  oceans  are  highly  time-dependent  regions.  Although 
observations  are  insufficient  to  determine  whether  the  region  near  the  western  bound- 
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800m  LOOPERS 


Figure  5.1:  A  collection  of  “looper”  float  trajectories  at  800  meters  in  the  Atlantic.  Notice  the 
development  of  eddy  structures  just  north  of  the  equator.  This  figure  is  reproduced  from  Richardson 
and  Schmitz  (1993). 

ary  in  the  tropics  is  m  nature  more  eddy  rich  than  at  mid-latitudes,  SOFAR  float 
measurements  in  the  Atlantic  {Richardson  and  Schmitz  (1993))  present  a  particularly 
compelling  example  that  suggests  this  possibility.  Shown  in  Figure  5.1,  is  a  set  of 
looper  trajectories  at  800  meters  depth  revealing  the  northward  propagation  of  ed¬ 
dies  along  the  boundary  m  this  region  and,  potentially,  their  generation  as  well. 

Close  examination  of  the  eddy  formation  process  in  the  numerical  model  reveals 
that  the  phenomenon  is  dependent  on  both  the  Reynolds  number  of  the  flow  and  also 
on  latitude.  The  eddies  begin  to  develop  for  nonlinear  flows  with  Re  >  Rcc,  where  the 
Reynolds  number  is  defined  by  global  basin  parameters  according  to  Equation  2.23. 
In  addition,  the  features  form  preferentially  near  the  equator  in  the  no-slip  case.  This 
behavior  is  illustrated  by  the  series  of  numerical  experiments  shown  in  Figure  5.2.  In 
all  three  cases,  the  source  is  located  at  the  northern  end  of  the  domain  and  the  sink 
m  the  southern,  and  the  basin  has  dimension,  A*  =  3  and  Ly  =  20  though  the  figures 
show  only  a  fraction  of  this  region  to  concentrate  on  the  unforced  central  region.  The 
system  is  spun  up  from  rest  to  a  steady  or  quasi-steady  state.  In  Figure  5.2a,  a  western 


boundary  current  carries  fluid  laminarly  from  the  source  region  across  the  equator  to 
the  sink  region.  Thus  the  flow  is  stable  at  this  Reynolds  number  {Re  =  25).  In  Figure 
5.2b,  the  Re  =  50  and  the  western  boundary  current  is  unstable,  and  eddies  develop. 
The  coincidence  between  their  formation  and  the  equator  is  striking  and  unmistakable. 
Finally,  Figure  5.2c  presents  a  run  at  Re  =  50  that  does  not  cross  the  equator.  Its 
central  latitude  is  at  t/  =  10,  and  the  flow  is  stable.  It  is  the  combination  of  a  high 
Reynolds  number  and  a  constraint  to  cross  the  equator  that  stimulate  the  observed 
eddy  field. 

Thinking  generally  about  the  stability  of  a  viscous  shallow-water  tropical  system 
suggests  three  mechanisms  that  could  be  important.  The  most  obvious  candidate 
is  horizontal  shear  instability,  which  has  an  extensive  history  dating  back  at  least  to 
Rayleigh  (1880)  (see  discussion  in  Drazin  and  Reid  (1981))  In  this  process,  the  velocity 
profile  of  the  mean  field  is  unstable  to  small  disturbances.  A  sufficient  condition  for 
the  flow  to  be  stable  is  the  absence  of  a  potential  vorticity  extrememum  in  the  flow 
profile  anywhere  in  the  domain.  As  a  result,  this  process  is  sometimes  referred  to 
as  an  inflection  point  instability.  Shear  instability  does  not  rely  on  the  presence  of 
friction.  To  the  contrary,  viscosity  that  is  present  in  the  system  acts  as  a  stabilizing 
force  which  can  completely  suppress  an  instability  if  sufficiently  strong. 

The  second  possibility  is  viscous  instability.  This  mechanism  has  been  suggested 
in  studies  of  the  mid- latitude  western  boundary  current  by  lerley  and  Young  (1991) 
and  Cessi  et  al.  (1993).  The  physics  of  this  instability  is  fundamentally  the  same  a.s 
above.  However,  the  basic  flow  (usually)  has  no  potential  vorticity  extremum  in  the 
domain  and  is  therefore  stable  to  perturbations  in  the  absence  of  friction.  Without 
viscous  forces,  the  Reynolds  stresses  {u'v')  vanish  owing  to  a  precise  phase  relationship 
between  the  different  velocity  components.  As  a  result  the  disturbance  can  not  extract 
momentum  from  the  mean  flow.  With  the  addition  of  viscosity,  this  phase  relationship 
is  broken  in  a  critical  layer  residing  in  the  high  shear  zone  adjacent  to  the  boundary, 
enabling  the  Reynolds  stresses  to  transfer  momentum  to  the  perturbation.  An  example 
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of  viscous  instability  is  the  plane  Pouseille  flow  which  has  a  parabolic  proflle.  With 
no  inflection  point,  the  flow  is  stable  to  shear  instability,  but  in  real  flows  at  very  high 
Reynolds  numbers  (Re=5780  (Kundu  (1990)),  the  system  undergoes  a  transition  to 
time-dependent  motion. 

A  third  possibility  is  that  the  flow  is  inertially  unstable,  a  process  more  directly 
related  to  the  rotation  of  the  system  than  horizontal  shear  instability.  A  familiar 
example  is  of  Taylor-Couette  flow  in  which  two  concentric  cylinders  with  fluid  between 
them  rotate  differentially.  It  can  be  shown  (e.g.,  Kundu  (1990)),  that  when  the  square 
of  the  circulation  of  fluid  decreases  with  radius,  the  system  becomes  unstable  and 
horizontal  rolls  develop  to  transfer  momentum  laterally.  Since  this  process  depends  on 
the  centrifugal  acceleration  of  the  fluid  parcels,  it  is  also  called  centrifugal  instability. 
In  geophysical  contexts  in  which  the  Coriolis  parameter  arises,  the  necessary  condition 
for  instability  in  a  highly  idealized  flow  is  that 

fq  <  0.  (5.1) 

When  the  potential  vorticity,  q,  of  the  flow  is  of  the  opposite  sign  as  the  Coriolis 
parameter,  as  happens  in  inertial  cross-equatorial  flow,  the  system  becomes  unsta¬ 
ble.  The  theory  applies  to  inviscid,  parallel  flow.  In  a  real  fluid,  friction  damps  the 
disturbance  as  in  shear  instability. 

In  nature,  other  instabilities  are  possible,  such  as  baroclinic  or  Kelvin- Helmholtz 
instability.  However,  these  depend  on  the  vertical  velocity  shear,  and  therefore  are 
absent  in  a  shallow-water  model.  Furthermore,  baroclinic  instability  is  not  likely 
for  mean  flows  with  horizontal  scale  small  with  respect  to  the  deformation  radius 
(Pedlosky  (1987a))  as  is  true  for  the  western  boundary  current  at  the  equator.  As  a 
result,  I  focus  only  on  the  three  mechanisms  above. 

Determining  the  physics  that  drives  the  equatorial  instability  is  not  a  pedantic 
quest.  Its  importance  rests  in  the  generality  of  the  eddy  formation  process  and  the 
dynamical  significance  attributed  to  the  eddy  field  in  the  previous  chapter.  If  eddies 
develop  only  as  a  result  of  the  details  of  the  frictional  parameterization,  as  suggested 
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by  the  possibility  of  a  viscous  instability,  then  their  role  in  transporting  vorticity 
to  the  western  boundary  current  is  also  quite  fragile,  potentially  just  an  artifact  of 
the  particular  model  chosen.  On  the  other  hand,  if  the  western  boundary  current  is 
found  to  be  unstable  to  the  more  robust  shear  instability,  then  the  conclusions  drawn 
previously  are  more  general  and  are  likely  to  apply  to  the  ocean  and  to  any  model 
that  produces  a  physically  realistic  western  boundary  current.  It  is  in  fact  difficult  to 
state  absolutely  which  process  is  fundamentally  responsible  if,  as  is  the  case,  all  the 
elements  for  each  instability  are  present.  However,  I  will  show  that  the  likely  process 
at  work  is  classical  shear  instability  of  the  western  boundary  current  that  is  triggered 

preferentially  at  the  equator  because  the  maximum  in  the  deformation  radius  occurs 
at  that  latitude. 

5.2  Formulation  of  the  Problem 

The  stability  of  the  tropical  western  boundary  current  is  assessed  using  the  standard 
techniques  of  linear  stability  analysis,  as  found  for  example  in  Pedlosky  (1987a)  or 
Drazm  and  Reid  (1981)  along  with  several  approximations.  In  particular,  the  western 
boundary  current  is  idealized  to  be  purely  meridional  and  latitude-independent  and 
the  effect  of  the  gradient  in  the  planetary  vorticity  is  neglected  in  the  perturbation. 
The  shallow-water  equations  in  dimensional  form  are  recast 


Ut  -f  u  •  Vu  —  fv  = 

—  +  S 

(5.2a) 

Ut  -f  u  •  Vu  -f  /u  = 

—  hy  +  +  T 

(5.2b) 

-I-  u  •  Vh  -f  hV  •  u  = 

(5.2c) 

where  additional  body  force  terms  c:(x,y)  and  T(x,y),  whose  utilility  will  become 
evident  below,  are  included  in  the  momentum  equations.  The  dependent  variables  are 
decomposed  into  time-mean  and  perturbation  quantities, 

u{x,y,t)  =  u{x,y)  +  u\x,y,t)  (5.3a) 
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v{x,y,t)  =  v{x,y)  +  v\x,y,t) 
h{x,y,t)  =  h{x,y)  +  h'{x,y,t). 


(5.3b) 

(5.3c) 


For  the  predominantly  meridional  shear  flow  of  a  western  boundary  current,  the  basic 
state  is  idealized  to  be  parallel  flow  that  is  only  a  function  of  longitude: 

u  =  0;  u  =  v(x);  h  =  h{x).  (5.4) 

Substituting  Equations  5.4  and  5.3a  —  5.3c  into  the  set  5.2a  —  5.2c  yields  a  new 
set  of  equations  in  terms  of  both  the  mean  and  perturbation  quantities.  Since,  linear 
stability  theory  models  the  initial  development  of  infinitesimal  disturbances  to  a  basic 
state,  the  equations  may  be  examined  at  different  orders  of  the  amplitude  of  the 
disturbance.  To  leading  order,  only  terms  independent  of  the  prime  variables  remain, 

fv  =  -g'hx  +  H  (5.5a) 

0  =  .<4//tJxi  +  T-  (5.5b) 

Equations  5.5a  and  5.5b  are  in  essence  a  redefinition  of  the  basic  state.  Body  force 
terms,  E  and  T,  remain  in  the  equations  at  this  order,  and  they  represent  the  artificial 
but  necessary  forces  to  balance  the  idealized  mean  flow  and  maintain  a  purely  merid¬ 
ional,  latitude-independent  basic  state.  Their  functional  form  depends  exactly  on  the 
chosen  form  for  the  mean  flow  used.  For  example,  for  a  desired  basic  state  velocity 
profile,  t/*(x)  and  ^ (x), 

H(x,y)  =  mA^)  +  g'K  (5-6a) 

T  =  -AhA-  (5-6b) 


Alternatively, 


dy 


=  /3u®. 


(5.7) 


Thus  the  forcing,  H,  satisfies  a  vorticity  balance  and  maintains  the  desired  meridional 
velocity  field,  and  T  balances  any  viscous  dissipation  of  that  flow.  The  vertical  mass 
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fltlx,  t«*,  is  absent  in  the  equations,  reflecting  the  fact  that  a  local  vertical  mass  flux 
can  not  drive  a  steady,  latitude-independent  basic  flow. 

Moving  to  the  next  order  retains  terms  linear  in  the  perturbations 

u't  -1-  vv! y  —  -f  h' X  —  Ah'^^u'  =  0  (5.8a) 

v't  +  Vxu’  +  vv'y  -f  fu'  -f  h'y  -  AhV^v'  =  0  (5.8b) 

^  t  +  hxU  -f  vh' y  -b  h  {u  X  +  v'y)  =  0.  (5.8c) 

The  time  evolution  of  the  disturbance  is  determined  by  this  set  of  equations  in  which 

the  perturbation  variables  are  coupled  both  to  one  another  and  to  the  basic  state.  The 
artificial  body  forces,  H  and  T,  that  drive  the  mean  flow  are  absent  in  the  perturbation 
equations.  Thus  their  inclusion  in  the  original  physics  has  no  effect  on  the  instability 
of  the  system  except  in  so  far  as  establishing  the  mean  field. 

The  coupled  set  of  partial  differential  equations,  5.8a-5.8c,  can  be  analytically 
simplified  by  making  a  stringent  assumption.  With  the  exception  of  the  Coriolis 
parameter,  /  =  /(y),  all  of  the  coefficients  of  the  dependent  variables  are  functions  of 
longitude  only  and  not  of  latitude  or  time.  If  &  local  /-plane  approximation  is  made 
and  the  Coriolis  parameter  is  assumed  to  have  a  constant  value,  /o,  then  a  modal 
decomposition  of  the  dependent  variables  is  possible  in  y  and  t.  Let 

(v',h')  =  (5.9a) 

u'  =  (5,9b) 

where  /  is  the  meridional  wavenumber  of  the  disturbance  and  c  represents  its  phase 
speed.  The  wavenumber-weighting  of  the  zonal  amplitude  function  in  Equation  5.9b 
merely  simplifies  the  ultimate  form  of  the  mathematical  problem  and  represents  no 
fundamental  difference  between  the  different  components.  The  method  outlined  here 
is  similar  to  that  of  Ni  (1996)  who  considered  a  zonal  flow  in  a  multi-layer,  inviscid 
shallow-water  model. 

The  /-plane  approximation  that  allows  this  decomposition  is  a  significant  limi¬ 
tation  of  the  present  analysis  as  the  variation  of  the  Coriolis  parameter  is  usually 
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essential  in  tropical  dynamical  studies.  A  full  treatment  of  the  problem  must  include 
both  the  meridional  variation  of  /  and  the  non~zonal  nature  of  the  basic  flow.  How¬ 
ever,  the  full  analysis  is  beyond  the  scope  of  this  work,  and  as  will  be  shown  below,  the 
simplification  to  a  local  /-plane  calculation  provides  surprisingly  good  agreement  with 
numerical  results  and  does  so  in  a  desirably  simple  and  understandable  framework. 

Using  relations  5.9a  and  5.9b,  the  system  5.8a--5.8c  reduces  to 

— (u  —  c)  U  —  UAh  {Uxx  ”  —  fV  +  Hx  =  0  (5.10a) 

(/  +  Vx)  U  +  (v  —  c)V  +  jAh  {Vxx  —  ^  ^  (5.10b) 

(hxU  +  hUxx)+hV  +  (v-c)H  =  0.  (5.10c) 

Equations  5.10a-5.10c  can  be  re-written  in  the  matrix  form 


Atp  =  ctj) 


(5.11) 


where 


and 


A  = 


UAh  {i£:  -  1 

)  i 

1  d 

12  dx 

f  +Vx 

V  +  UAh  (p  5^2  “ 

1  1 

hx  +  hA 

h 

V 

V’  = 


u 

V 

H 


(5.12) 


(5.13) 


is  the  state  vector.  Equation  5.11  is  a  standard  eigenvalue  problem  in  which  there  exists 
a  corresponding  eigenvector,  for  each  eigenvalue,  c*.  From  the  modal  representation 
of  the  perturbation  amplitudes,  stable  modes  have  phase  speeds  either  real  or  complex 
with  negative  imaginary  component.  Any  eigenmode  that  has  a  positive  imaginary 
component  grows  exponentially  in  time,  at  least  until  the  perturbations  reach  finite 
amplitude. 

Equation  5.11  gives  the  stability  properties  for  meridional  shear  flow,  of  infinite 
zonal  extent.  To  apply  this  formulation  to  the  realistic  shear  flow  of  a  western  bound¬ 
ary  current  requires  a  finite  domain,  x  =  [O,®^;],  where  the  distance  to  the  eastern 
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boundary  of  the  basin,  xe,  is  expected  to  be  large  compared  to  the  scale  of  the  shear 
flow  itself.  Boundary  conditions  at  the  solid  meridional  boundaries  are  no  flux  and  no 
slip, 

=  0  (5.14a) 

=  0.  (5.14b) 

The  full  linear  stability  analysis  even  for  this  somewhat  idealized  flow  as  expressed 
by  the  matrix,  A,  reveals  a  complicated  dependence  of  the  state  vector  on  the  basic 
state  and  on  the  latitude,  through  the  magnitude  of  the  Coriolis  parameter.  It  is  not 
immediately  clear  how  the  stability  of  the  system  will  change  with  latitude  or  with  the 
mean  flow  itself.  To  gain  further  insight  from  the  equations  directly  requires  additional 
simplifications  that  are  discussed  in  the  next  section. 


5.3  Inviscid  Shear  Instability 


Although  Equation  5.11  can  be  discretized  and  analyzed  as  written  for  a  variety  of 
mean  fields,  it  is  useful  to  further  reduce  the  problem  first  so  as  to  isolate  the  different 
instability  mechanisms.  Setting  the  viscosity  coefficient,  to  zero  eliminates  viscous 
instability  from  the  system.  The  matrix,  A  then  reduces  to 


A  = 


f  +  Vx  V  1 


(5.15) 


[hx  +  h^  h  u  J 

and  the  order  of  the  system  is  reduced.  The  new  eigenvalue  problem  is  defined  by. 


AV’  =  cip-  (5.16) 

Only  two  boundary  conditions  are  needed,  ajad  the  obvious  choice  is  that  of  no  normal 
flow.  Equation  5.14a.  The  no-slip  condition.  Equation  5.14b  seems  unphysical  if  the 
disturbance  can  not  detect  any  dissipation.  Thus  friction  is  required  to  establish  the 
mean  flow  profile,  but  its  effect  on  the  perturbation  is  neglected. 
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The  stability  criterion  for  inviscid  zonal  flows  on  the  j8-plane  in  the  shallow-water 
system  has  been  examined  by  Ripa  (1983).  The  case  with  =  0  corresponds  to  the 
present  system.  Stability  is  insured  if  (i) 

[a  —  t7(a;)]  ^  >  0  (^-1^) 

ax 

and  (a), 

[a  —  v{x)f  <  g'h{x)  (5.18) 

for  all  X  and  any  a.  For  the  standard  Munk  boundary  current  as  well  as  inertio- viscous 
profiles,  condition  5.17  is  generally  violated  assuming  a  in  the  range  of  the  velocity. 

The  following  section  briefly  describes  the  numerics  involved  in  the  solution  of 
the  eigenvalue  problem.  Section  5.3.2  determines  and  examines  the  stability  of  5.16 
using  the  inviscid  formulation,  A,  and  shows  that  the  Munk  boundary  layer  is  in  fact 
unstable  to  small  perturbations,  and  that  the  shear  instability  growth  rate  is  largest 
near  the  equator.  Section  5.3.3  makes  two  additional  approximations  that  allow  a 
still  simpler  mathematical  and  physical  interpretation  of  the  phenomenon  that  more 
clearly  illustrates  the  tropical  enhancement  of  the  instability. 


5.3.1  Numerics 

All  eigenvalue  equations  analyzed  are  discretized  on  a  stretched  grid  in  x  with  stag¬ 
gered  variables  in  a  manner  similar  to  that  described  in  Section  A. 2  for  the  time- 
dependent  model,  except  that  for  this  western  boundary  problem,  the  stretching  func¬ 
tion  is  a  simple  exponential: 

x(i)  =  d  (e'i^  —  ij  (5.19) 


where 


e*  -  1’ 


(5.20) 


N  is  the  number  of  grid  points,  and  s  gives  the  degree  of  stretching.  A  value  of 
s  =  3  WcLs  chosen  to  adequately  enhance  the  resolution  in  the  rapidly  varying  western 
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boundary  region  without  introducing  excessive  deformation  of  the  grid.  The  results 
are  slightly  dependent  on  the  number  of  grid  points.  All  calculations  were  determined 
using  a  value  of  at  least  N  =  100,  though  even  N  =  Z2  appears  adequate  to  reveal 
the  general  behavior  of  the  instability. 

Generally,  the  matrix.  A,  is  determined  numerically,  given  the  basic  state  and  the 
latitude,  and  the  eigenvalues,  c,  and  eigenvectors,  'ip,  are  determined  using  the  eigen¬ 
value  solver,  eig.m,  which  calls  industry  standard  EISPACK  routines,  and  is  included 
in  the  basic  distribution  of  Matlab.  Checking  the  solutions  using  a  different  solver, 
sptarn.m,  which  applies  a  different  algorithm  to  find  selected  eigenvalues  produces 
virtually  indistinguishable  results.  The  numerical  method  for  solution  was  checked 
against  and  agreed  well  with  the  asymptotic  stability  analysis  of  Lipps  (1963)  for  a 
wall-free  Bickley  jet. 


5.3.2  Instability  of  the  Munk  Boundary  Layer 


It  is  desirable  to  examine  first  a  simple  basic  flow  in  which  the  fields  are  truly  latitude- 
independent.  Any  asymmetry  that  develops  in  the  instability  characteristics  can  then 
be  attributed  to  the  shallow-water  system  and  not  to  changes  in  the  potential  vorticity 
structure  of  the  basic  state.  I  will  examine  later  in  Section  5.3.4  the  instability  of 
a  more  realistic,  latitude-dependent  flow  field.  The  analytic  solution  to  the  Munk 
boundary  layer  provides  a  convenient  basic  velocity  profile  to  begin  the  investigation: 


v(x)  =  Dsin 


(5.21) 


where 


D  = 


250 


\/3^. 


M 


(5.22) 


Rather  than  use  the  geostrophically  balanced  height  field,  which  varies  with  latitude, 
I  set  the  height  of  the  layer  uniform: 


h  =  l. 


(5.23) 
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a)  Munk  velocity  profile 


0.029  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9 


X 

Figure  5.3:  a)  Velocity  (in  units  of  c  =  y/g'H)  and  b)  potential  vorticity  gradient  (in  for 
the  Munk  boundary  layer  with  6m  =  0.1  and  Sq  =  —0.025.  Only  shown  to  a:  =  1,  though  xe  =  2 
in  calculations  below.  Dot-d2ished  curve  marks  the  intersection,  v  —  0.059,  which  equals  the  phase 
speed  of  the  fastest  unstable  mode  with  1  =  6. 

The  quantitative  adjustment  for  a  geostrophically  balanced  height  field  for  the  weak 
velocities  examined  is  small.  The  Munk  layer  is  determined  using  =  0.1,  and  for 
a  transport  of  So  =  0.025,  consistent  with  the  numerical  model  runs  in  Figure  5.2. 
Equation  5.16  is  solved  for  the  profile,  5.21,  a  Coriolis  parameter,  /  =  1,  and  a  reduced 
gravity,  g*  =  1.  The  basic  state  velocity  profile  is  shown  in  Figure  5.3  along  with  the 
zonal  gradient  of  the  potential  vorticity.  The  eastern  boundary  extends  20  6m j  though 
only  the  westernmost  10  6m  is  shown  in  the  figure.  The  potential  vorticity  gradient 
passes  through  zero  at  multiple  longitudes  so  the  inflection  point  criterion  is  satisfied 
by  this  basic  state. 

Note  that  the  scales  of  the  basic  states  used  in  this  chapter  are  identical  to  those 
of  the  rest  of  the  thesis.  Thus  length  scales  are  given  in  terms  of  deformation  radii 
(i.e.,  wavenumbers  in  terms  of  inverse  deformation  radii),  growth  rates  are  scaled  by 
PLd,  and  phase  speeds  are  relative  to  the  gravity  wave  speed,  c  =  \/g'H. 
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The  growth  rate,  Im(/c)  =  w.,  of  all  eigenmodes  of  the  system  are  shown  in  Figure 
5.4a  as  a  function  of  the  wavenumber,  1.  The  number  of  modes  equals  the  dimension, 
ZN,  of  the  system.  Most  modes  are  neutral,  with  zero  imaginary  component.  They 
represent  disturbances  that  are  purely  wavelike  that  pass  through  the  system  neither 
growing  nor  decaying.  Positive  and  negative  indicate  growing  and  decaying  modes, 
respectively,  and  the  symmetry  about  Ici  —  0  indicates  that  for  each  growing  mode 
there  exists  a  corresponding  decaying  mode.  This  feature  is  to  be  expected  from  the 
form  of  the  matrix  5.15  in  which  all  elements  are  purely  real.  The  state  vector,  can 
be  replaced  by  its  complex  conjugate,  and  the  equation  remains  the  same.  Thus  there 
is  a  symmetry  between  eigenvalues  and  their  complex  conjugates,  and  therefore  also 
between  growth  and  decay  rates.  Notice  that  this  feature  changes  with  the  inclusion 
of  friction,  as  in  matrix.  A,  Equation  5.12. 

In  Figure  5.4a  a  few  sets  of  unstable  modes  can  be  traced  as  a  function  of  wavenum¬ 
ber.  In  particular,  one  pair  of  modes  has  the  largest  amplitude  for  all  /,  and  they  peak 
at  a  value  of  approximately  I  =  3.25.  This  series  and  the  additional  sets  of  modes  that 
appear  within  its  envelope  are  the  sets  of  modes  that  I  consider  to  result  from  shear 
instability,  and  they  are  quite  insensitive  to  the  exact  discretization  and  grid  resolu¬ 
tion.  Notice  that  at  high  wavenumber  there  exist  an  increasing  number  of  unstable 
modes  of  small  amplitude.  These  are  more  fragile  modes  that  depend  more  closely  on 
the  discretization  of  the  system.  While  all  unstable  modes  grow,  it  is  the  one  with  the 
largest  growth  rate  that  dominates  the  instability.  Therefore  these  numerically  sen¬ 
sitive  modes  can  be  neglected  when  considering  also  the  more  unstable  shear  modes. 
Furthermore,  their  amplitude  tends  to  decrease  with  increasing  resolution.  Within  the 
approximations  of  the  analysis.  Figure  5.4a  indicates  that  the  Munk  layer  is  unstable 
to  inviscid,  growing  perturbations. 

The  phase  speed  of  the  fastest  growing  mode  is  shown  as  a  function  of  wavenumber 
in  Figure  5.4b  The  value  varies  between  a  value  of  =  —0.98  at  /  =  0.01  and 
Cr  =  —.056  at  /  =  6,  both  in  the  range  of  v. 
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The  amplitudes  of  the  eigenvectors  (e.g.,  y/UU*)  associated  with  the  fastest  grow¬ 
ing  modes  are  shown  in  Figure  5.5.  The  zonal  velocity  and  height  fields  are  smooth 
functions  of  longitude,  with  a  noticeable  maximum  between  x  —  0.3  and  x  =  0.5, 
depending  on  the  wavenumber.  At  high  wavenumbers,  the  eigenfunctions  show  the 
rise  of  a  second  peak  near  x  =  0.1.  The  V  component  also  has  an  interesting  two  peak 
structure,  with  a  third,  particularly  sharp  feature  developing  for  high  wavenumbers. 
The  eigenfunctions  are  very  smooth  for  long  wave  solutions.  The  sharpest  peak  at 
high  wavenumber  is  associated  with  the  critical  location  where  Munk  velocity  profile 
equals  the  phase  speed  of  the  fastest  growing  mode, 

v{x)  -  cl.  (5.24) 

For  a  jet  profile  with  the  phase  speed  in  the  range  of  the  velocity  field,  this  relation 
is  satisfied  at  least  at  two  longitudes.  From  Figure  5.4  it  is  clear  that  v  =  0.059 
corresponds  to  the  phase  speed  of  the  fastest  growing  mode  at  /  =  6.  This  line  is 
overlayed  as  a  dot-dashed  curve  in  Figure  5.3a.  Equality  5.24  applies  at  the  intersection 
of  the  two  curves.  The  first  intersection  occurs  at  x  =  0.275,  which  corresponds  closely 
to  the  peak  in  the  eigenvector  amplitude  at  x  =  0.029.  The  second  intersection  occurs 
at  X  =  .262,  where  the  potential  vorticity  gradient,  is  nearly  zero.  Thus  the 
asymmetry  of  the  Munk  boundary  layer  appears  to  play  a  role  in  the  development 
of  the  strong  peak  in  the  eigenvector  curve.  It  occurs  approximately  at  the  location 
determined  by  Equation  5.24,  but  the  one  which  does  not  satisfy 

^  ^  0.  (5.25) 

ax 

Identical  analyses  for  a  range  of  latitudes  provide  the  growth  rate  as  a  function 
of  wavenumber  and  latitude.  Letting  f  —  fiy  with  =  1  provides  the  growth  rate 
explicitly  as  a  function  of  latitude.  This  field  is  contoured  in  Figure  5.6a  for  the 
same  Munk  layer  used  above.  Solid  lines  separate  intervals  of  a  hundredth  and  two 
dot  dashed  contours  at  values  of  uJi  =  0.058  and  0.0584  are  used  to  aid  the  eye.  A 
growth  rate  maximum  =  0.0584  resides  clearly  in  the  immediate  vicinity  of  the 
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Figure  5.6:  As  functions  of  meridional  wavenumber  and  latitude,  the  a)  growth  rate  and  b)  phase 
speed  of  the  fastest  growing  mode  for  the  Munk  boundary  layer  in  Figure  5.3.  The  dot-dashed  curves 
in  (a)  mark  w,  =  0.058  and  0.0584. 

equator,  at  a  wavenumber  I  =  3.2  and  at  a  latitude  of  y  =  —0.5.  However,  this  growth 
rate  maximum  is  not  located  at  y  =  0.  Since  the  basic  state  is  independent  of  y,  it 
appears  that  the  shallow-water  system  is  aware  of  the  direction  of  the  flow  through 
the  absolute  vorticity,  /  +  Ux,  in  matrix,  A. 

The  phase  speed  of  the  fastest  growing  mode  is  contoured  in  Figure  5.6b  also  as 
a  function  of  wavenumber  and  latitude.  For  a  given  wavenumber,  the  phase  speed 
of  the  most  unstable  disturbance  decreases  with  latitude.  The  significance  of  this 
structure  along  with  a  more  physical  interpretation  of  the  equatorial  enhancement  of 
the  instability  is  further  explored  in  the  context  of  a  still  simpler  model  in  the  next 
section. 
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5.3.3  Geostrophic  Perturbations 

The  previous  section  shows  that  the  Munk  boundary  layer  is  unstable  to  inviscid  per¬ 
turbations,  and  that  the  instability  is  greatest  near  the  equator.  However,  the  compli¬ 
cated  form  of  matrix.  A,  prevents  a  direct  explanation  for  the  tropical  enhancement. 
Two  assumptions  that  follow  examples  of  Stem  (1961)  and  Lipps  (1963)  significantly 
simplify  the  system.  First,  the  mean  flow  is  assumed  only  weakly  nonlinear  so  that 

(5.26) 

This  approximation  is  clearly  invalid  very  close  to  the  equator,  where  /  =  0,  but  the 
flow  can  be  specified  to  be  arbitrarily  weak  so  that  relationship  5.26  holds  for  all  other 
latitudes. 

Second,  the  perturbations  must  be  geostrophically  balanced, 

-/«'  =  (5.27a) 

fu'  —  —g'h'y.  (5.27b) 

This  relation  assumes  that  the  Rossby  number  of  the  perturbation  is  small,  which 
is  valid  again  arbitrarily  close  to  the  equator  for  truly  infinitesimal  motion.  It  also 
assumes  the  time  rate  of  change  of  the  perturbation  be  small  relative  to  the  Coriolis 
term, 

Ic 

y  <  1.  (5.28) 

Since  the  growth  rate  of  a  disturbance  in  linear  theory  scales  with  the  velocity,  this 
relation  can  again  be  satisfied  arbitrarily  close  to  the  equator  for  sufficiently  weak 
flows. 

The  geostrophic  approximation  couples  the  different  fields  of  the  disturbance  such 
that  a  single  equation  results,  expressed  in  terms  of  a  single  variable.  Adopting  the 
above  approximations  and  scaling  the  horizontal  dimensions  by  the  length  scale,  L, 
yields  for  the  height  anomaly,  H, 

(u  -  c)  +  f)h]-  -Fv)H  =  0,  (5.29) 
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where 


(5.30a) 

(5.30b) 


is  the  rotational  Fronde  number.  At  the  boundaries,  the  height  of  the  disturbance  is 
set  to  zero: 


H  =  0. 


(5.31) 


Equation  5.29  describes  the  initial  development  of  the  instability  for  weak  mean 
flows  until  the  Rossby  number  of  the  disturbance,  becomes  order  1.  For  F  =  0, 
Equation  5.29  reduces  to  the  familiar  Rayleigh  equation  or  the  equation  for  standard 
barotropic  instability.  The  additional  terms  proportional  to  F  reflect  the  influence 
of  the  horizontal  divergence  in  the  velocity  field  that  is  included  in  the  shallow-water 
model. 

Separating  the  amplitude,  H ^  and  the  phase  speed,  c,  into  real  and  imaginary  parts 
and  integrating  over  the  domain  gives  the  usual  necessary  criterion  for  the  development 
of  an  instability, 

G\H\^  , 

^  /  7= - =  0,  (5.32) 

Jo  |u  —  c|  ^  ^ 

where 

G  =  v^^-  Fv  (5.33) 

is  the  mass- weighted  potential  vorticity  gradient.  Equation  5.32  shows  that  for  growing 
modes  to  exist  (i.e.,  if  c^  >  0),  there  must  be  a  location  within  the  domain  where  G 
changes  sign. 

The  solution  to  Equation  5.29  subject  to  boundary  conditions  5.31  gives  remarkably 
similar  results  to  that  of  the  full  system.  The  contour  plot  of  growth  rate  as  a  function 
of  wavenumber  and  latitude  is  shown  in  Figure  5.7a.  Solid  contours  are  uniform,  every 
hundredth,  and  the  dot-dashed  contours  correspond  to  values,  =  0.058  and  0.0584 
as  in  Figure  5.6a.  The  maximum  of  0.0588  is  located  at  a  wavenumber  of  Z  =  3.26 
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Phase  speed 


Figure  5.7:  As  functions  of  meridional  wavenumber  and  latitude,  the  a)  growth  rate  and  b) 

phase  speed  for  the  fastest  growing  mode  as  in  Figure  5.6,  but  assuming  geostrophic  perturbations 
as  in  Equation  5.29.  The  wavenumber  corresponding  to  the  growth  rate  maximum  as  a  function  of 
latitude  is  shown  as  the  dashed  curve  in  both  plots.  The  dotted  curve  is  the  wavenumber,  Imj  defined 
in  Equation  5.36  as  a  function  of  latitude. 

and  y  =  0.  Note  that  there  is  no  asymmetry  in  this  modified  instability  problem  as 
Equation  5.29  can  not  distinguish  between  the  northern  and  southern  hemispheres  if 
the  basic  state  is  longitude  independent. 

The  effect  of  the  horizontal  divergence  of  the  velocity  field  can  now  be  interpreted 
in  Equation  5.29  more  simply  than  through  the  full  shallow-water  system.  Introducing 
the  effective  wavenumber,  k  such  that 

^  =  l^  +  F,  (5.34) 

Equation  5.29  becomes 

-  k^H  -  :^H  =  0.  (5.35) 

V  —  c 

Consider  first  the  ca.se  where  G  and  v  are  independent  of  latitude.  In  this  case, 
Equation  5.35  has  exactly  the  form  of  the  Rayleigh  equation  in  terms  of  the  effective 
wavenumber,  k.  A  generic  growth  rate  curve  for  an  unstable  mean  flow  is  shown  in 
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Figure  5.8:  Generic  growth  rate  curve  for  a  wavenumber,  /c. 


Figure  5.8.  The  growth  rate  has  a  maximum  at  a  value  of  which  is  independent  of 
latitude.  As  the  Coriolis  parameter  increases  with  latitude,  the  deformation  radius  Ld 
decreases,  and  therefore  F  increases.  The  actual  meridional  wavenumber  associated 
with  Kjj^  is 

(5.36) 

which  decreases  rapidly  with  latitude.  However,  the  actual  growth  rate  depends  on  I 
and  not  «,  so  the  true  growth  rate  maximum  does  not  follow  a  relation  like  equation 
5.36,  but  rather  a  more  modest  decrease  with  latitude.  Both  and  the  actual  growth 
rate  maximum,  Im,  are  shown  in  Figure  5.7  as  functions  of  latitude. 

This  trend  can  be  expressed  alternatively.  For  a  fixed  actual  wavenumber,  Zq, 
increasing  F  with  latitude  shifts  k  to  a  larger  value.  If  Iq  <  Im,  this  effect  can  induce 
larger  growth  rates  at  that  wavenumber,  initially.  However,  this  increase  ends  once 
is  reached.  Further  increases  in  latitude  reduce  the  growth  rate  of  the  instability  until 
ultimately  the  high  wavenumber  cutoff,  /ci,  is  passed,  and  the  flow  becomes  stable. 

Next  examine  the  effect  of  divergence  through  the  last  term  on  the  left  hand  side 
of  5.35.  From  the  stability  criterion.  Equation  5.32,  it  is  clear  that  if  G  is  of  one  sign 
then  Ci  =  0.  If  the  flow  has  an  inflection  point  and  F  is  small,  then  this  requirement 
is  satisfied.  However,  if  the  flow  is  of  one  sign,  as  for  example  with  a  fully  inertial 
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boundary  current,  then  at  some  critical  latitude,  F(yc}v  >  for  all  x,  and  the  flow 
is  stable.  In  this  case,  the  effect  of  divergence  is  also  stabilizing.  For  flows  that  vanish 
at  least  once  in  the  domain,  such  as  the  Munk  boundary  layer,  the  flow  can  satisfy 
the  necessary  criterion  for  instability  at  all  latitudes,  but  the  flow  may  become  stable 
at  high  latitudes  according  to  the  previous  argument. 

Thus  the  flow  becomes  more  stable  as  the  ratio  of  the  length  scale  of  the  mean 
flow  to  the  deformation  radius  increases.  A  physical  interpretation  is  obtained  by 
considering  the  instability  as  the  result  of  interacting  vorticity  anomalies. 

Any  curvature  in  the  velocity  profile  acts  as  a  vorticity  gradient  upon  which  Rossby 
waves  propagate.  The  phase  speed  of  a  wave  in  a  meridional  mean  flow  is  given  in 
terms  of  the  deformation  radius  by 


(5.37) 


where  the  effective  vorticity  gradient  is 


=  -V, 


(5.38) 


A  schematic  of  the  western  boundary  current  is  presented  in  Figure  5.9,  in  which  it 
is  clear  that  the  vorticity  gradients  change  sign  on  either  side  of  the  inflection  point, 
P.  Relative  to  the  velocity,  tJp,  of  the  mean  flow  at  P,  waves  travel  southward  to  its 
east  and  northward  to  its  west,  but  are  advected  by  the  relative  mean  flow.  Therefore 
it  is  possible  to  have  waves  that  are  stationary  with  respect  to  vp.  These  phase- 
locked  modes  interact,  advecting  anomalous  vorticity  and  encouraging  the  growth  of 
neighboring  disturbances. 

As  the  deformation  radius  increases,  the  phase  speed  of  a  given  wave  increzises.  To 
maintain  stationarity  with  respect  to  the  mean  flow,  the  wavenumber  must  decrease. 
However,  decreasing  I  increases  the  separation  between  anomalies,  and  they  interact 
less  strongly.  Actual  growth  is  determined  by  the  balance  of  phase-locking  and  mutual 
reinforcement.  If  only  stationarity  were  important,  the  wavenumber,  Im,  associated 
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Figure  5.9:  A  western  boundary  current  profile  with  an  inflection  point,  P.  Note  that  the  potential 
vorticity  gradient  due  to  the  current  changes  sign  at  P  which  acts  as  the  local  "north”  pole. 


with  the  maximum  in  the  growth  rate  would  follow  the  dotted  curve  in  Figure  5.7b 
which  shows  Note  its  close  association  with  the  lines  of  constant  phase.  In  fact, 
decreases  much  less  rapidly  with  latitude  to  maintain  the  necessary  proximity  of 
anomalies. 

The  physics  of  the  instability  can  also  be  interpreted  from  the  perspective  of  ener- 
getics.  In  the  context  of  this  simplified  treatment  that  considers  only  geostrophically 
balanced  perturbations,  it  is  possible  to  derive  the  energy  equation, 


(5.39) 


Here,  the  overbar  symbol  with  a  y  superscript,  implies  an  averaging  in  the  t/-direction. 
Zonally  integrated,  the  time  evolution  of  total  energy  of  the  disturbance  depends 
on  the  ability  of  the  disturbance  to  extract  energy  from  the  mean  flow  through  the 
Reynolds  stress,  Tl^H^ . 

The  total  energy  is  divided  between  potential,  the  last  term  on  the  left  hand  side 
of  Equation  5.39,  and  kinetic,  the  first  two  terms,  and  their  ratio  is  0(F).  At  mid- 
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latitudes,  F  is  large,  and  most  of  the  energy  of  the  disturbance  is 

partitioned  into 

potential  form  (P).  In  contrast,  near  the  equator  the  fraction  of  total  energy  that  is 
kinetic  [K)  reaches  its  maximum.  However,  from  the  right  hand  side  of  Equation  5.39, 
it  is  clear  that  the  growth  of  the  instability  is  proportional  to  the  Reynolds  stress,  and 

the  Schwarz  inequality  allows 

,  ,  (H^f  4-  (Hyf 

{H.Hy)  <  ^  2  "• 

(5.40) 

Seeking  solutions. 

H^H{x,y)e^\ 

(6.41) 

it  can  be  shown  that 

jdx{c(K  +  P)<K(v,)^J, 

(6.42) 

or 

(6.43) 

Thus  the  growth  rate  of  the  disturbance  depends  on  the  ratio  of  the  kinetic  energy  to 
the  total  energy.  When  the  partitioning  favors  potential  energy,  this  ratio  is  small  and 
the  disturbance  grows  more  slowly  than  when  more  of  the  energy  is  in  kinetic  form.  At 
the  equator,  where  the  system  is  most  barotropic  and  the  kinetic  energy  dominates, 
the  disturbance  extracts  energy  from  the  mean  flow  most  efficiently.  At  mid-latitudes, 
where  much  of  the  energy  is  stored  in  layer  height  variations,  the  Reynolds  stress  is 
small,  and  the  instability  is  more  subdued. 

It  is  important  to  note  that  the  original  work  of  Lipps  (1963)  illustrated  the  sta¬ 
bilizing  effect  of  divergence  on  a  zonal  mid-latitude  flow,  and  this  work  was  extended 
to  the  equator  by  Philander  (1976). 

5.3.4  Instability  of  an  Inert io- Viscous  Current 

The  analysis  above  reveals  that  the  Munk  boundary  current  is  inviscidly  unstable. 
However,  western  boundary  currents  found  in  nature  and  in  realistic  numerical  models 


139 


likely  do  not  satisfy  a  linear  vorticity  balance.  To  examine  the  stability  of  a  more 
realistic  mean  state,  the  output  of  the  time-dependent  numerical  model  is  used  to 
provide  a  basic  state.  The  model  output  is  shown  in  Figure  5.2a,  which  is  spun  up  just 
below  threshold  with  Re  =  25.  Note  that  this  profile  technically  violates  the  parallel 
flow  assumption,  5.4,  because  it  varies  slightly  with  latitude.  However,  the  variations 
with  longitude  are  much  larger  than  those  in  y,  and  this  detail  is  ignored.  Velocity  and 
height  profiles  are  taken  directly  from  the  numerical  model  output  and  interpolated 
onto  the  relevant  numerical  grid  for  the  instability  calculations  using  a  cubic  spline 
interpolant.  This  section  applies  the  full,  inviscid  shallow- water  formulation  using  A 
as  opposed  to  the  simplified  model  of  the  last  section. 

Figure  5.10a  presents  the  two-dimensional  plot  of  growth  rate  versus  latitude  and 
meridional  wavenumber.  The  highest  wavenumber  data  have  been  cropped  to  omit 
noise  that  appears  there  from  the  numerical  modes  that  slightly  dominate  there.  Com¬ 
parison  with  Figure  5.6  show  that  the  two  profiles  have  distinctly  similar  instability 
characteristics,  despite  the  addition  of  a  mildly  latitude-dependent  basic  state  and 
a  slightly  more  nonlinear  profile.  The  maximum  growth  rate,  =  0.0574,  is 

slightly  less  than  that  for  the  Munk  profile,  and  it  occurs  at  I  =  3.2  and  y  =  —.5 
identical  to  the  other  calculation.  The  phase  contours  m  Figure  5.10b  are  also  very 
similar  to  the  Munk  curves,  though  slightly  shifted  in  magnitude  apparently  by  slight 
differences  in  profile. 

It  is  surprising  that  the  behavior  of  the  latitude-independent  and  latitude-varying 
profiles  are  so  similar.  The  potential  vorticity  gradient  of  a  single  layer  model  takes 
the  form 


dq 

dx 


(5.44) 


If  the  height  field  is  in  geostrophic  balance  with  the  velocity,  then  the  denominator  in 
fhe  rightmost  term  above  might  distinguish  the  stability  properties  of  one  profile  versus 
the  other  and  induce  its  own  hemispheric  asymmetry.  However,  the  term  is  multiplied 
by  p,  so  its  effect  is  small  near  the  equator  where  the  dynamics  are  dominated  by  the 
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Mode!  output:  a)  Growth  rate 
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b)  Phase  speed 
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Figure  5.10:  As  functions  of  meridional  wavenumber  and  latitude,  the  a)  growth  rate  and  b) 
phase  speed  of  the  fastest  growing  mode  as  in  Figure  5.6  except  that  the  basic  state  used  is  from 
the  numerical  model  experiment  shown  in  Figure  5.2a.  The  dot-dashed  curve  in  (a)  demarcates 
Wi  =  0.057. 


first  term.  Thus  the  instability  properties  are  quite  robust  to  meridional  variations  in 
the  basic  state. 

5.4  The  Influence  of  Viscosity 

The  previous  calculations  applied  to  inviscid  perturbations.  However,  viscosity  is  an 
essential  part  of  the  dynamics  of  the  mean  flow,  and  therefore  it  should  be  expected  to 
influence  the  dynamics  of  the  perturbation.  Indeed  the  fact  that  the  system  is  stable 
for  some  Reynolds  numbers  shows  that  its  effect  is  not  negligible.  The  eigenvalue 
problem  including  dissipation  along  with  the  extra,  no-slip  boundary  condition  is 
solved  similarly  to  the  above  inviscid  problems.  The  stability  of  the  basic  state  is 
examined  for  a  range  of  viscosities. 

The  eigenvalues  for  numerical  model  output  data  at  y  =  1  determined  using  a 
viscosity.  Ah  =  2.5  x  {Re  =  100)  are  shown  in  Figure  5.11.  One  feature  that 
distinguishes  these  results  from  those  in  Figure  5.4  is  that  the  eigenvalues  no  longer 
come  in  complex  conjugate  pairs.  The  introduction  of  viscosity  changes  the  symmetry 
of  the  problem,  and  causes  most  modes  to  decay  in  time.  In  fact,  only  those  modes 
near  the  imaginary  axis  are  shown  in  the  Figure  so  as  to  illustrate  the  more  interesting 
single  set  of  growing  modes.  These  correspond  directly  to  those  fastest  growing  modes 
in  the  inviscid  problem,  now  slightly  damped  by  the  introduction  of  a  small  amount  of 
friction,  but  not  sufficiently  to  damp  the  instability.  One  change  is  that  the  maximum 
m  the  growth  rate  curve  is  shifted  out  to  higher  wavenumber,  /  =  4.  In  addition, 
the  marginal  points  that  bound  the  region  of  stability  from  instability  are  shifted 
somewhat  from  the  inviscid  problem.  The  phase  speeds  associated  with  that  set  are 
shown  in  Figure  5.11b,  and  have  a  similar  character  to  those  of  the  inviscid  problem. 

The  amplitudes  of  the  eigenfunctions  of  the  unstable  modes  are  shown  in  Fig¬ 
ure  5.12  and  again,  their  similarity  to  the  inviscid  case  is  also  striking.  One  distinguish¬ 
ing  feature  of  the  viscous  modes  is  found  in  the  amplitudes  of  the  V  eigenfunctions.  As 
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Mode!  output.  Viscous  perturtatioo  (Re=100):  a)  Growth  rate 


Figure  5.11:  As  functions  of  meridional  wavenumber,  the  a)  growth  rate  and  b)  phase  speed  of  the 
fastest  growing  mode  for  the  viscous  stability  problem  using  the  model  output  profile.  Calculations 
assume  Re  =  100  to  compare  directly  with  the  inviscid  calculations  of  the  Munk  boundary  layer. 
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in  the  previous  problem,  a  peak  appears  at  small  x  in  the  functions  corresponding  to 
large  wavenumber.  However,  the  peak  is  broader  than  in  the  inviscid  case  suggesting 
that  viscosity  plays  an  important  role  in  the  dynamics  of  the  critical  layer. 

These  calculations  apply  a  large  Reynolds  number  explicitly  to  compare  directly 
the  inviscid  calculations  of  the  Munk  boundary  layer  to  an  approximately  inviscid 
calculation  with  a  small  amount  of  viscosity.  The  close  agreement  between  modes 
and  eigenfunctions  suggests  strongly  that  the  basic  mechanism  of  the  instability  is 
identical  to  the  inviscid  problem,  namely  shear  instability. 

Calculations  at  several  latitudes  show  similar  bullseye  patterns  of  growth  enhanced 
near  the  equator,  but  with  reduced  amplitude  depending  on  the  magnitude  of  the 
viscosity  employed.  For  example,  using  a  viscosity  of  8.333  x  10“^  which  corresponds 
to  Re  =  30  gives  a  maximum  growth  rate,  =  0.014,  centered  at  I  =  3.8  and 

y  =  —1-2.  Increasing  the  dissipation  appears  to  drive  the  instability  further  south  of 
the  equator,  and  lengthen  the  time  needed  for  its  development.  The  shift  to  a  higher 
wavenumber  is  found  consistently  at  all  latitudes. 

This  highly  idealized  viscous  theory  indicates  that  the  numerical  model  should  be 
unstable  for  Re  =  25.  In  fact,  numerical  integrations  indicate  that  the  instability 
begins  near  Re  =  32.  This  small  discrepancy  may  result  from  the  approximations  of 
the  instability  calculation  itself  (e.g.,  neglect  of  the  variation  in  /).  The  numerical 
model  may  in  fact  be  unstable  for  Re  —  25.  The  time-scale  for  growth  of  the  instability 
at  this  Reynolds  number  is  0(5000)  time  units,  which  is  a  very  long  integration  that 
I  have  not  examined.  At  Re  =  32  the  time-scale  is  0(300)  time  units,  which  is  the 
time-scale  over  which  they  are  observed  to  develop.  Thus  the  numerical  model  may 
be  unstable  at  smaller  Reynolds  numbers  than  Re  =  32,  but  so  slightly  that  the 
time-dependent  motion  escapes  detection. 

A  second  measure  of  the  accuracy  of  the  instability  calculation  is  to  compare  the 
length  scales  of  the  disturbance  predicted  by  the  theory  to  that  found  in  the  numerical 
experiment.  After  reaching  an  approximate  steady  state  at  Re  —  31.25,  the  viscosity 
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Figure  5.12:  Eigenfunction  amplitudes  of  a  subset,  /  =  1  to  i  =  6  in  intervals  of  0.5,  of  the  fastest 
growing  modes  in  Figure  5.11:  (a)  U;  (b)  V;  (c)  H.  Notice  the  broaxl  peak  in  the  V  eigenfunction 
near  x  =  0.05 


Figure  5.13:  The  height  field  anomaly  after  instability  begins.  The  difference  in  height  fields 
between  steady  solution  at  Re  =  31.25  to  Re  =  32.25  after  300  time  units. 

of  the  numerical  model  was  changed  slightly  to  Re  =  32.25,  and  integrated  further 
in  time.  Figure  5.13  presents  the  difference  in  the  height  fields  after  300  time  units 
have  elapsed.  Height  anomalies  appear  just  south  of  the  equator  as  expected  from 
the  theory,  and  develop  further  to  the  south.  There  is  some  variability  being  swept 
eastward  from  the  western  boundary  near  the  southern  boundary  of  the  plot,  where 
the  boundary  current  veers  eastward  as  part  of  the  southern  hemisphere  sink.  Approx¬ 
imately  one  full  wave  is  visible  in  the  figure,  and  it  has  a  wavelength  of  approximately 
1.5  deformation  radii.  This  wavelength  agrees  well  with  that  predicted  by  the  theory 
of  A  =  Thus  even  with  the  extensive  idealizations  of  the  theory,  it  predicts  the 
behavior  of  the  disturbances  found  in  this  numerical  model,  in  location  for  growth,  in 
wavelength,  and  in  the  critical  Reynolds  number. 

The  numerical  example  described  above  contains  eddies  that  are  approximately  of 
the  deformation  scale,  consistent  with  the  linear  theory  and  the  basic  state  applied. 
To  demonstrate  that  the  eddy  dimension  is  in  fact  determined  by  the  width  of  the 
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aty=-3(-),y=0(-.).y=3(~) 


Figure  5.14:  a)  Velocity  and  b)  height  fields  as  a  function  of  x  at  three  latitudes  in  the  basin 
(y  =  3  (“),  y  =  0  and  y  —  —Z  (-)).  This  profile  for  a  second  experiment  has  a  zonal  scale  that  is 
roughly  half  that  shown  in  Figure  5.2a. 

boundary  current,  and  not  the  deformation  radius,  I  have  examined  a  second  numerical 
experiment  in  which  the  viscosity  has  been  reduced  by  a  factor  of  eight.  Setting 
Ah  =  1.25  X  10  a  steady  solution  results  from  using  H  =  1  and  5^0  =  0.003.  For 
these  parameters.  Re  =  24  is  just  below  the  critical  value  at  which  time-dependent 
motion  begins,  which  indicates  that  Rec  is  not  a  universal  number,  but  varies  with 
model  parameters,  such  as 

Inserting  the  numerical  solution  at  time  to  =  1000  (shown  at  three  latitudes  in 
Figure  5.14)  into  the  viscous  stability  model  as  its  basic  state  gives  growth  rates  as 
a  function  of  the  viscosity  that  affects  the  perturbation.  For  a  Reynolds  number  of 
24,  the  growth  rate  is  nearly  zero.  Choosing  a  value  for  the  viscosity.  Ah  =  1.071  x 
10”^,  reveals  that  the  boundary  current  is  unstable  to  infinitesimal  perturbations 
at  Re  =  28.  The  growth  rates  and  phase  speeds  for  this  experiment  are  shown 
in  Figure  5.15  as  functions  of  meridional  wavenumber  and  latitude.  Note  that  the 
growth  rate  maximum  occurs  at  a  meridional  wavenumber  of  /  =  7.6,  approximately 
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Model  output,  Viscous  perturbation  (Re=28):  a)  Growth  rate 


Figure  5.15:  As  functions  of  meridional  wavenumber  and  latitude,  a)  growth  rates  and  b)  phase 
speeds  of  the  most  unstable  modes  for  the  basic  state  shown  in  Figure  5.14,  using  a  viscosity  of 
1.071  X  10  (iZe  28  for  the  perturbations).  The  dashed  contour  corresponds  to  a  value  of 
=  10  PLd  and  divides  the  domain  into  regions  of  growth  and  regions  of  stability. 

twice  the  value  in  the  former  experiment  (see  again  Figure  5.11)  and  consistent  with 
the  smaller  boundary  current  here.  The  bullseye  for  this  run  is  also  somewhat  more 
extended  meridionally  from  the  previous  examples,  showing  that  as  the  boundary 
current  narrows  with  respect  to  the  deformation  radius  and  the  system  becomes  more 
barotropic,  the  tropical  enhancement  to  the  instability  is  less  pronounced. 

To  compare  the  growth  rates  predicted  by  the  theory  to  those  observed  in  the 
model,  I  carry  out  further  integration  of  this  numerical  experiment  at  the  lower  vis¬ 
cosity  {Ah  =  1.071  X  10“^).  This  exercise  generates  cyclonic  structures  in  the  vicinity 
and  south  of  the  equator  much  as  in  Figure  5.2b,  but  of  a  smaller  scale.  Showing  the 
zonal  velocity  anomaly,  du  =  u(to  +  At)  -  u{to),  at  intervals  of  At  =  400,  Figure  5.16 
presents  the  time-evolution  of  the  instability.  The  contour  interval  is  the  same  in  each 
plot  (du  =  0.002c).  After  the  first  time  interval,  the  western  boundary  current  has 
superposed  on  it  a  series  of  very  weak  structures  of  alternating  sign.  After  At  =  800 


148 


(/3Lu)"\  well-defined  eddies  of  alternating  sign  emerge  between  y  =  -1  and  y  =  -5. 
As  time  continues,  the  variability  grows,  in  strength  and  in  meridional  extent.  After 
1600  time  units,  the  eddies  exist  north  of  the  equator,  and  the  strength  of  those  near 
y  =  —3  have  peaked  in  strength. 

Although  the  negative  anomalies  (dashed  in  the  figure)  tend  to  be  slightly  larger  in 
magnitude  than  the  positive,  there  is  very  good  symmetry  between  positive  and  neg¬ 
ative  features.  This  remains  true  for  the  meridional  velocity  and  height  fields,  though 
the  vorticity  field  shows  a  more  pronounced  asymmetry  favoring  positive  anomalies. 
This  symmetry  supports  the  modal  decomposition  of  the  u,  v,  and  h  fields  in  Equa¬ 
tions  5.9a  and  5.9b.  The  time-scale  for  the  eddy  development  in  the  model  is  800 
to  1000  {PLd)~^,  slightly  less  than  the  value  predicted  by  the  viscous  theory  (1300 
()0i/£))“^),  and  somewhat  longer  than  that  predicted  by  the  inviscid  calculation  (250 
{I3Ld)~^,  not  shown).  The  fully  developed  anomalies  in  Figure  5.16d  have  a  char¬ 
acteristic  wavelength  of  between  |  and  1  deformation  radius,  depending  on  the  lat¬ 
itude  at  which  the  estimate  is  made.  The  value  predicted  by  the  stability  analysis 
is  A  =  =  0.8.  The  inviscid  calculation  predicts  a  wavelength  of  approximately  1 

deformation  radius.  The  theory  reasonably  predicts  the  time  and  space  scales  of  the 
eddies,  their  location  for  development,  and  the  critical  Reynolds  number  at  which  they 
appear. 

The  addition  of  viscosity  to  the  system  has  adjusted  the  instability  properties 
slightly,  shifting  the  wavenumber  at  which  the  maximum  growth  rate  is  found  and 
reducing  the  amplitude  of  the  growth.  However,  it  has  not  produced  a  qualitative 
change  in  the  development  of  the  system.  Viscosity  acts  as  a  stabilizing  influence,  but 
the  physics  of  the  instability  is  unaltered.  These  features  suggest  that  the  instability 
mechanism  responsible  for  the  growth  of  the  eddies  in  the  time-dependent  model  is 
classical  shear  instability  and  not  viscous  instability. 
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Figure  5.16:  Zor 

of  the  eddy  field  afi 
from  24  to  28. 


5.5  Inertial  Instability 


Before  concluding  that  shear  instability  is  solely  responsible  for  the  development  of 
the  eddies,  it  is  worthwhile  to  examine  briefly  the  possibility  of  inertial  instability. 
As  mentioned  in  the  introduction  to  this  chapter,  this  instability  arises  in  geophysical 
flows  when  the  local  absolute  vorticity  is  of  the  opposite  sign  to  the  local  planetary 
vorticity,  as  is  likely  to  occur  in  cross-equatorial  flow.  This  section  derives  the  necessary 
stability  criterion  for  inertial  instability  in  the  shallow  water  system,  but  shows  that 
it  is  a  not  likely  to  be  acting  on  the  western  boundary  current  of  the  model. 

To  derive  the  instability  criterion,  return  again  to  the  inviscid  system,  described 
by  Equation  5.16.  In  order  to  isolate  inertial  instability  from  shear  instability  it  is 
helpful  to  examine  a  basic  flow  structure  that  has  no  potential  vorticity  extremum. 
The  analysis  in  this  section  thus  considers  flows  that  differ  substantially  from  the 
western  boundary  currents  discussed  above  to  illustrate  the  basic  mechanism  of  inertial 
instability.  In  particular,  I  examine  the  stability  of  parabolic  flows  which  vary  very 
slowly  in  longitude. 

The  WKB  expansion  (e.g..  Bender  and  Orszag  (1978))  begins  with  the  introduction 
of  a  small  parameter,  e,  which  defines  a  slow  variable,  X  —  ex,  over  which  changes 
in  the  character  of  the  eigenfunction  solution  occurs.  Eigenfunctions  are  expressed  in 
ternis  of  x  and  X  as 

U(x)  =  U{X)e^  (5.45) 

and  similarly  in  V  and  H.  The  phase,  k,  is  also  assumed  a  slowly  varying  function 
and  can  be  expressed  as  a  series  in  powers  of  the  small  parameter,  e: 

K^)  =  E 

n=l 

The  amplitude  and  phase  of  the  perturbation  can  vary  on  the  scale  x,  but  the  rate  of 
change  of  the  phase  varies  slowly.  Then 

U.  =  -j:e-^U  +  eU^e^  (5.47a) 
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(5.47b) 


,  (ik\  __  ^  /  V 

To  leading  order,  Equation  5.16  reduces  to 


V  —  c 


f  +  Va 
hj.  “f"  h/i- 


V  —  c 


(5.48) 


A  non-trivial  solution  requires  that  the  matrix  in  Equation  5.48  have  a  non-zero  de¬ 
terminant.  Then, 

fT-  A  [/'-  ^2  r  .  /(/  +  ^x)  i  dki  (-  .dfci-'\l  I  ( r-r  _  .dki-\ 

c)  1^(7;  c)  +  _  __  q.  ^ 

(5.49) 

which  gives  the  general  stability  properties  for  the  shallow-water  system  under  the 
assumptions  of  the  WKB  approximation  (i.e.,  that  the  rate  of  change  of  the  phase  of 
the  disturbance  is  small). 

The  conditions  that  imply  c  complex  are  not  obvious  from  the  form  of  Equa¬ 
tion  5.49.  As  a  result,  I  consider  the  limit  of  — >  0,  which  requires  that  the 
amplitudes  of  the  eigenfunctions  vary  according  to  the  slow  variable,  X.  Physically 
this  assumption  should  correspond  to  the  disturbance  of  a  very  slowly  varying  mean 
field.  If  the  disturbance  senses  a  homogeneous  basic  state,  then  the  eigenfunction 
must  be  approximately  independent  of  location.  Assuming  also  that  the  mean  height 
and  velocity  fields  are  geostrophically  balanced,  and  examining  only  that  particular 
location  where  u  =  0,  then  three  eigenvalues  are  derived. 


—  di\  h -\- 


f(f  +  Vx) 


(5.50a) 

(5.50b) 


The  associated  frequencies  are 


a;  =  0,  ±^Ph  +  f{f  +  v^). 


(5.51) 


For  any  /(/  -f  v^)  <  0,  there  is  some  wavenumber  for  which  the  radical  in  Equation 
5.51  is  negative,  and  u  is  complex.  The  maximum  growth  rate  occurs  for  1  =  0,  which 
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distinguishes  this  process  from  shear  instability  whose  maximum  occurs  for  non-zero 
1.  Thus  inertial  instability  in  the  two-layer  fluid  generates  long- wave  disturbances. 
Requiring  oj  real  in  Equation  5.51  gives 

f{f  +  vs)>0  (5.52) 


which  is  exactly  the  stability  criterion  derived  by  Ripa  (1983)  for  a  zonal  flow  on  the 
/5-plane  assuming  purely  longitude-independent  motions  (/  =  0  in  this  case). 

Relation  5.52  is  independent  of  conditions  5.17  and  5.18  and  simply  provides  the 
necessary  condition  that  ensures  inertial  stability.  Absolute  stability  is  guaranteed  by 
the  criteria  5.17  and  5.18.  A  flow  may  be  found  to  violate  the  condition  for  inertial 
instability,  5.52,  but  satisfy  the  conditions  for  the  shallow-water  system,  and  therefore 
be  stable. 

To  see  whether  the  criterion  5.52  has  a,ny  practical  use,  I  have  applied  the  inviscid 
stability  equation,  5.16,  to  a  flow  that  is  expected  to  be  inertially  unstable.  The  mean 
flow  has  a  parabolic  profile  of  amplitude  P, 


7 


(x  -  Xo)^'\ 
T2  )  ’ 


(5.53) 


where  xq  =  f  is  the  center  of  the  profile  in  the  center  of  the  channel  of  width,  L,  and 


7  =  4  (5.54) 

ensures  the  velocity  goes  to  zero  at  the  boundaries.  Without  an  inflection  point,  this 
profile  isolates  the  experiment  from  shear  instability.  The  basic  state  height  field  is  in 
geostrophic  balance  with  the  velocity. 

Setting  P  =  2  violates  criterion  5.18  and  using  /o  =  0.01  ensures  that  condition 
5.52  is  not  satisfied.  Calculations  from  this  flow,  with  no  stretching  and  N  =  100,  show 
no  evidence  of  a  long- wave  instability.  Similar  runs,  with  different  parameters  appear 
to  be  stable  also.  In  fact,  I  have  been  unable  to  find  any  instability  in  a  parabolic 
profile  on  an  /-plane  except  for  those  in  which  the  potential  vorticity  gradient  changes 
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sign,  which  can  occur  if  relation  5.54  is  relaxed,  and  the  velocity  is  allowed  to  change 
sign.  ^ 

While  the  stability  criteria  5.17,  5.18  and  5.52  are  suggestive,  they  provide  no 
insurance  that  the  flows  are  in  fact  unstable.  Inertial  instability  appears  to  be  very 
fragile,  depending  sensitively  on  the  details  of  the  jet.  Returning  to  the  subject  of  the 
instability  of  the  western  boundary  current,  it  appears  unlikely  that  inertial  instability 
is  responsible  for  the  development  of  the  eddies.  Not  only  is  the  physical  mechanism 
sensitive,  but  it  also  predicts  long-wave  features,  in  contrast  to  the  roughly  deformation 
radius  scale  eddies  that  develop  in  the  model. 


5.6  Summary  and  Discussion 

Linear  stability  analyses  of  parallel  flow  in  the  /-plane  shallow- water  system  show  that 
the  western  boundary  current  is  unstable  to  infinitesimal  perturbations,  and  that  this 
instability  is  greatest  in  the  vicinity  of  the  equator.  Three  possible  driving  mecha¬ 
nisms  are  shear,  viscous,  and  inertial  instability.  Inviscid  shear  instability  calculations 
agree  well  with  numerical  model  results,  predicting  both  the  correct  length  scale  for 
the  disturbances  and  the  latitude  where  they  develop.  Viscosity,  when  added  to  the 
theory,  acts  primarily  to  damp  the  growth  of  the  instability  and  modifies  the  eigen¬ 
functions  of  the  growing  modes  only  quantitatively,  and  in  the  vicinity  of  a  critical 
layer.  These  facts  lead  me  to  conclude  that  viscous  instability  is  not  the  dominant 
process.  The  viscous  calculations  do  predict  the  critical  Reynolds  number  at  which  the 
western  boundary  current  becomes  unstable  in  reasonable  agreement  with  the  time- 
dependent  experiments.  An  analytic  treatment  of  the  shallow-water  system  suggests 
that  it  supports  inertial  instability  under  very  restrictive  conditions,  but  numerical 
calculations  of  realistic  flow  structures  (though  not  resembling  the  western  boundary 

^Interestingly,  the  parabolic  profile  with  vanishing  boundary  velocity  does  appear  to  be  unstable 
also  for  a  zonal,  e^terly  jet  on  the  /3-plane  centered  on  the  equator.  This  profile  is  also  distinguished 
in  that  the  meridional  potential  vorticity  gradient  changes  sign. 


154 


current),  show  that  this  type  of  instability  is  not  robust.  It  seems  unlikely  that  this 
process  plays  an  important  role  in  the  instability  of  the  western  boundary  current. 

These  results  are  consistent  with  mid-latitude  model  results  of  lerely  and  Young  (1991). 
That  study  includes  the  meridional  variation  of  the  Coriolis  parameter,  but  they  also 
examine  the  case  of  /3  =  0.  However,  while  they  determine  that  the  central  process 
is  viscous  instability,  I  conclude  that  the  boundary  current  suffers  a  shear  instability. 
This  statement  implies  that  the  results  of  this  study  are  more  generally  applicable. 
Any  boundary  current  whose  velocity  scale  decreases  or  simply  levels  off  to  a  constant 
value  near  the  boundary,  such  as  the  mean  velocity  in  Figure  1.3,  satisfies  this  criterion 
and  therefore  is  likely  unstable  provided  a  sufficiently  high  Reynolds  number. 

Consistent  with  numerical  model  results,  the  tropics  exhibit  preferred  eddy  growth 
in  both  the  viscous  and  inviscid  calculations.  The  latitude  dependence  results  from 
the  increase  of  the  internal  deformation  radius  with  decreasing  latitude  which  modifies 
the  phase  speeds  of  the  waves  that  travel  on  the  vorticity  gradient  of  the  mean  field. 

It  is  ironic  that  it  is  the  baroclinic  nature  of  the  shallow-water  model  that  allows  the 
equator  to  distinguish  itself  from  all  other  latitudes,  and  yet  it  is  the  fact  the  local 
dynamics  at  the  equator  are  most  barotropic  that  allows  for  the  enhanced  instability 
there. 

Several  approximations,  such  as  the  neglect  of  the  meridional  variation  in  the  plan¬ 
etary  vorticity  and  mean  flow,  have  aided  the  above  analysis  and  contributed  to  the 
simple  physical  interpretation.  It  is  possible  in  the  context  of  the  simple  model  of 
Section  5.3.3  to  relax  one  simplification  and  consider  the  effect  of  yS.  Carrying  out 
the  same  analysis  with  the  vorticity  equation  and  maintaining  the  divergence  in  the 
velocity  field  yields  a  similar  eigenvalue  relation  to  Equation  5.11.  Calculations  based 
on  this  model  show  similar  features  to  those  in  the  /-plane  model,  with  an  equa- 
torially  enhanced  growth  rate,  and  comparable  phase  speeds  of  the  unstable  modes 
(Figure  5.17).  Interestingly,  the  growth  rate  maximum  is  somewhat  larger  than  with¬ 
out  (compare  =  0.085  to  0.058).  As  in  the  study  of  lerley  and  Young  (1991), 


155 


Munk  layer,  Geostrophic  perturbations,  beta=1:  a)  Growth  rate 


Figure  5.17:  As  functions  of  meridional  wavenumber  and  latitude,  the  a)  growth  rate  and  b) 
phase  speed  as  in  Figure  5.7  except  that  /3  =  1  is  included. 

basin  modes  develop  at  the  small  wavenumbers,  a  result  of  the  meridional  vorticity 
gradient,  but  these  have  small  growth  rates  compared  to  the  modes  that  decay  with 

X. 

The  effect  of  a  vorticity  gradient  raises  interesting  possibilities  when  considering 
more  realistic  oceanic  applications.  For  example,  western  boundary  currents  are  sel¬ 
dom  found  along  purely  meridional,  perfectly  vertical  walls.  Borrowing  the  results 
from  studies  of  barotropic  instability  in  zonal  flows  suggests  that  a  topographic  /3, 
which  is  likely  to  be  relevant  for  deep  boundary  currents,  stabilizes  the  instability 
at  mid-latitudes  but  not  at  the  equator.  Thus  topographic  influences  would  further 
enhance  the  tropical  nature  of  the  instability.  Also,  the  effect  of  a  sloping  coast  breaks 
the  symmetry  of  the  instability  process,  with  westward  currents  being  more  unstable. 

The  agreement  between  the  simplified  analysis  and  the  numerical  results  is  some¬ 
what  remarkable,  given  the  strong  variation  of  the  planetary  vorticity  in  the  region. 
Ultimately,  a  more  formal  approach  to  the  problem  that  includes  the  meridional  vari- 
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ation  of  /,  the  non-zonal  nature  of  the  flow  and  potentially  several  other  features  is 
warranted.  In  addition,  an  analysis  of  the  free-slip  system  would  further  support  or 
negate  the  results  of  the  present  theory.  However,  the  physics  associated  with  the 
meridional  variation  in  the  deformation  radius  is  quite  robust,  and  I  anticipate  that 
shear  instability  will  remain  enhanced  at  the  equator. 
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Chapter  6 


Conclusion 


6.1  Summary  and  Discussion  of  Thesis 

This  thesis  examines  the  dynamics  of  nonlinear  cross-equatorial  flow  in  the  context  of 
idealized  western  boundary  currents  in  the  deep  ocean.  The  dynamical  problem  centers 
on  the  transformation  of  potential  vorticity  that  must  occur  for  fluid  parcels  to  merge 
smoothly  with  mid-latitude  circulations  in  the  opposite  hemisphere.  I  have  found 
that  the  system  behavior  is  qualitatively  divided  between  two  regimes.  The  linear, 
or  weakly  nonlinear  regime  is  characterized  by  steady  flow,  with  potential  vorticity 
modification  carried  out  by  frictional  processes  in  a  laminar  boundary  layer.  This 
regime  has  been  well-studied  by  previous  theory.  At  a  particular  ratio  of  nonlinearity 
to  dissipation  defined  by  a  critical  Reynolds  number,  the  system  passes  a  transition 
to  non-steady  flow  in  which  eddies  form  in  the  western  boundary  current  and  are 
advected  southward  across  the  equator  where  ultimately  they  dissipate.  This  thesis 
identifies  the  dynamical  role  of  the  eddy  field  as  the  necessary  conduit  of  vorticity  that 
allows  cross-equatorial  flow.  As  the  Reynolds  number  increases  beyond  the  critical 
value,  a  turbulent  boundary  layer  emerges  that  transports  vorticity  from  the  inertial 
portion  of  the  boundary  current  to  the  frictional  sub-layer  that  exists  adjacent  to 
the  boundary  where  potential  vorticity  modification  is  straightforward.  In  addition. 
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this  thesis  determines  that  the  cause  of  the  eddy  field  is  an  instability  of  the  western 
boundary  current  that  is  preferentially  triggered  at  the  equator  where  the  internal 
deformation  radius  is  at  its  maximum.  The  mechanism  for  the  instability  is  suggested 
to  be  classical  shear  instability  damped  by  the  presence  of  viscosity. 

Through  a  scaling  of  ErtePs  potential  vorticity  from  available  observations,  Chap¬ 
ter  2  determines  that  an  appropriate  model  for  the  study  of  cross-equatorial  flow  is 
the  shallow-water  model.  This  system  has  the  necessary  stretching  term  in  the  vor¬ 
ticity  equation  to  distinguish  the  equator  from  other  latitudes  and  provides  a  simple 
framework  for  the  study  of  cross-equatorial  flow.  Dissipation  is  included  as  viscous 
terms  in  the  momentum  equations. 

To  model  deep  flow  in  the  ocean,  the  classical  model  of  Stommel  and  Arons  (1960a) 
is  adopted  in  which  fluid  is  extracted  from  the  model  uniformly  throughout  the  domain, 
which  is  rectangular  in  geometry  and  straddles  the  equator.  To  balance  the  mass  sink 
in  the  interior  of  the  ocean  is  a  source  region  located  in  the  north-western  corner  of 
the  domain.  This  model  necessarily  constrains  fluid  in  the  western  boundary  current 
to  cross  the  equator. 

Chapter  3  introduces  and  describes  the  details  of  the  particular  numerical  ex¬ 
periments.  The  dynamics  of  the  western  boundary  current  are  determined  by  the 
Reynolds  number.  For  small  iZe,  the  system  is  stable  and  asymptotes  to  a  steady 
state.  The  dynamical  description  of  this  regime  basically  follows  the  original  Stommel 
and  Arons  (1960a)  solution  with  a  Munk  layer  appended  at  the  western  boundary.  For 
Reynolds  numbers  just  above  a  threshold  value,  the  boundary  current  becomes  unsta¬ 
ble  and  eddies  develop  just  north  of  the  equator.  The  eddies  propagate  steadily  across 
the  equator  in  a  periodic  or  quasi-periodic  fashion.  Only  cyclonic  eddies  emerge  at  this 
level  of  Reynolds  number,  reflecting  the  high  potential  vorticity  of  the  source  region. 
The  time-dependent  region  is  confined  to  approximately  one  eddy  length  scale  from 
the  western  boundary  current.  As  the  system  becomes  still  more  nonlinear,  vorticity 
anomalies  of  both  signs  develop,  and  the  domain  of  time-dependent  activity  extends 
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quite  far  into  the  interior  several  deformation  radii  north  and  south  of  the  equator. 
There  also  appears  substantial  coupling  of  energy  into  the  equatorial  waveguide  at 
this  level  of  forcing. 

Despite  the  intense  eddy  activity,  the  time-mean  circulation  remains  qualitatively 
unchanged  from  the  linear  picture.  There  exists  an  ever  intensifying  western  boundary 
current  that  transports  required  volume  of  fluid  across  the  equator  as  determined  by 
the  distribution  of  the  upwelling.  Regardless  of  the  Reynolds  number,  the  system  is 
able  to  respond  to  carry  fluid  across  and  change  its  potential  vorticity. 

Vorticity  analyses  of  the  time-mean  circulation  illustrate  the  dynamical  role  of 
the  eddy  field  in  modifying  potential  vorticity  and  are  presented  in  Chapter  4.  The 
method  examines  the  flux-conservative  form  of  the  time-mean,  absolute  vorticity  (i.e., 
the  mass-weighted  potential  vorticity)  equation  in  which  the  vorticity  flux  vectors 
are  subdivided  into  three  parts,  the  flux  of  mean  absolute  vorticity  by  the  mean 
velocity  field,  an  eddy  flux  term,  and  the  flux  due  to  viscosity.  Since  the  vorticity 
flux  vector  is  non-divergent,  it  is  straightforward  to  follow  vorticity  within  the  domain 
and  observe  that  a  substantial  fraction  of  vorticity  advected  southward  in  the  western 
boundary  current  exits  the  basin  through  the  western  boundary.  This  process  is 
explained  in  the  context  of  the  linear  and  weakly  nonlinear  limits.  Examination  of 
the  different  flux  terms  in  various  regions  of  the  domain  highlight  the  essential  physics 
operating  in  those  regions  and  illustrate  the  transition  from  one  mechanism  of  vorticity 
transfer  to  another.  In  the  laminar  experiments,  mean  absolute  vorticity  is  advected 
southward  by  the  western  boundary  current  and  is  transported  westward  by  the  viscous 
flux.  Thus  potential  vorticity  modification  is  simply  a  result  of  the  western  boundary 
current  being  largely  a  frictional  boundary  current.  At  higher  Reynolds  numbers 
where  inertial  effects  become  important,  the  viscous  boundary  layer  recedes  relative 
to  the  full  boundary  layer  and  the  westward  flux  is  carried  through  by  the  eddy  term. 
There  appears  a  turbulent  boundary  layer  which  grows  with  Reynolds  number  and 
transports  vorticity  from  the  relatively  inertial  portion  of  the  boundary  current  to  the 
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narrow  frictional  sub-layer. 

Numerical  experiments  carried  out  at  very  large  Reynolds  numbers,  Re  =  1000 
indicate  that  the  magnitude  of  the  eddy  flux  term  is  not  bounded  by  vorticity  input 
into  the  system.  As  the  system  becomes  increasingly  nonlinear,  the  eddies  must  com¬ 
pensate  for  an  increasingly  narrow  frictional  boundary  layer  and  for  an  eastward  flux 
of  vorticity  carried  by  the  mean  field.  This  behavior  contrasts  with  the  convective 
problem  in  which  the  eddy  flux  is  bounded  by  the  buoyancy  forcing  applied  to  the 
system. 

One  feature  that  emerges  in  the  numerical  experiments  at  high  Reynolds  number 
is  a  recirculation  in  the  vicinity  of  the  equator.  This  recirculation  is  not  associated 
with  trapped  fluid  of  constant  potential  vorticity  as  are  highly  nonlinear  recirculations 
^id-latitude  models.  Instead  it  appears  to  be  a  region  of  enhanced  eddy  formation 
or  perhaps  a  region  of  retarded  eddy  propagation.  This  eddy-driven  feature  has  the 
dynamical  role  of  contributing  to  the  necessary  westward  eddy  flux  of  vorticity  in 
this  important  region  just  north  of  the  equator.  Thus  both  the  mid-latitude  and 
equatorial  recirculations  are  essential  dynamical  elements  in  the  vorticity  balance,  but 
their  character  is  quite  different. 

The  role  of  the  eddy  field  is  quite  insensitive  to  details  of  the  model.  Numerical 
experiments  with  different  upwelling  distributions  and  different  boundary  conditions 
present  a  fundamentally  similar  impression.  Regardless  of  the  volume  of  fluid  that  is 
required  to  penetrate  into  the  opposite  hemisphere,  the  eddies  are  able  to  transport 
the  necessary  potential  vorticity  to  permit  the  flow.  The  boundary  conditions  merely 
establish  approximately  the  latitude  where  the  zonally  integrated  vorticity  flux  must 
vanish.  For  a  no-slip  condition,  it  is  the  equator,  and  for  free-slip,  it  is  south  of  the 
equator,  but  not  indefinitely  south. 

Finally,  Chapter  5  investigates  the  cause  of  the  eddy  held  itself.  Since  the  dynam¬ 
ics  of  the  cross-equatorial  flow  rest  fundamentally  on  the  existence  of  an  eddy  field  to 
transport  vorticity  westward  to  the  boundary,  it  is  desirable  to  understand  the  mech- 
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anism  that  creates  them.  A  linear  stability  analysis  of  an  idealized  western  boundary 
current  indicates  that  the  flow  is  unstable  to  inviscid  shear  instability.  The  boundary 
velocity  has  an  inflection  point  in  its  profile  and  therefore  does  not  satisfy  the  suffi¬ 
cient  conditions  for  stability.  Adding  viscosity  to  the  analysis  damps  the  growth  of  the 
instability,  but  does  not  qualitatively  alter  its  features.  The  meridional  wavenumber 
associated  with  the  maximum  growth  rate  changes  slightly,  but  the  eigenfunctions  of 
the  instability  are  nearly  identical. 

An  interesting  feature  that  emerges  from  the  numerical  model  as  well  as  from 
observations  and  the  analysis  is  that  the  instability  of  the  western  boundary  current 
is  tropically  enhanced.  The  maximum  in  the  growth  rate  occurs  in  the  vicinity  of 
the  equator.  A  slight  asymmetry  between  hemispheres  results  from  the  fact  that 
the  absolute  vorticity  of  the  western  boundary  current  is  not  symmetric  about  the 
equator.  The  physical  explanation  for  this  equatorial  enhancement  centers  on  the 
increased  scale  of  the  equatorial  deformation  radius.  The  meridional  variation  in  the 
deformation  radius  shifts  the  phase  speeds  of  waves  that  propagate  on  the  vorticity 
gradient  of  the  mean  flow  such  that  the  adjacent  disturbances  at  the  equator  have  the 
largest  wavenumber  and  can  mutually  reinforce  one  another  most  eff^ectively. 

It  is  important  to  reiterate  the  motivation  for  this  work.  Observations  indicate 
that  fluid  travels  between  hemispheres,  necessarily  across  the  equator  in  deep  western 
boundary  currents.  Although  observations  are  insufficient  to  determine  the  mecha¬ 
nism  for  vorticity  transformation,  it  is  certain  that  at  some  point  along  the  current 
path,  the  vorticity  is  modified.  If  the  deep  western  boundary  currents  are  friction- 
ally  dominated,  then  the  classical  studies  of  the  abyssal  circulation  are  sufficient  to 
explain  the  dynamics.  However,  scaling  arguments  from  observations  indicate  that 
the  boundary  currents  are  at  least  partially  inertial,  in  which  case  the  present  study 
applies  directly. 

The  model  predicts  the  formation  of  time-dependent  features  at  a  critical  Reynolds 
number.  Variability  is  observed  in  many  measurements  of  the  tropical  oceans,  at  all 
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depths.  A  dominant  time-scale  for  the  variability  is  from  10  to  70  days.  The  model  also 
finds  variability  on  these  time-scales.  It  is  the  scale  of  the  eddy  features  that  naturally 
emerge  in  the  model.  However,  these  model  features  can  be  tuned  by  parameter 
settings,  and  without  clear  understanding  of  the  lateral  viscosity  of  the  ocean,  it  is 
difficult  to  make  direct  comparisons  other  than  to  say  that  the  sets  of  features  are 
analogous  and  probably  generated  identically. 

Based  on  the  above  calculations  of  critical  Reynolds  numbers,  it  is  possible  to 
estimate  the  critical  lateral  viscosity  at  which  the  western  boundary  goes  unstable. 
Using  the  definition  Re  =  let  Rcc  =  30.  Then  for  a  deep  western  boundary 

current  that  is  approximately  1000  m  deep  and  that  transports  20  Sv  across  the 
equator,  the  critical  viscosity  at  which  variabity  would  set  in  would  be  Ah  =  670m^ s~U 

6.2  Extension  to  Full  Three-Dimensional  Fluid 

It  is  tempting  to  explore  the  process  of  cross-equatorial  flow  using  a  more  realistic, 
fully  three-dimensional  ocean  model.  Of  course,  such  an  experiment  would  retain 
substantial  dynamical  constraints  that  are  omitted  using  the  shallow-water  model. 
However,  I  believe  that  the  basic  description  described  above  would  remain  essentially 
unaltered  by  such  a  study. 

Fundamentally,  the  vorticity  flux  constraint  remains  identical  whether  describing 
a  two-  or  three-dimensional  fluid.  The  formalism  of  the  J- vector  analysis  clarifies  this 
point.  The  impermeability  theorem  states  that  mass- weighted  potential  vorticity,  pq^ 
can  not  be  transported  across  isentropic  surfaces  (i.e.,  surfaces  of  constant  potential 
density)  {Haynes  and  McIntyre  (1987)).  Thus  the  three-dimensional  fluid  is  conve¬ 
niently  decomposed  into  a  series  of  layers,  each  one  analogous  to  the  shallow-water 
system  described  in  this  thesis. 

The  form  of  the  vorticity  flux  vector  takes  a  slightly  different  form  for  a  three- 
dimensional  fluid.  In  the  notation  of  Marshall  and  Nursur  (1992),  the  potential  vor- 
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ticity  equation  can  be  written 


^(ot)  +  vj  =  o 


(6.1) 


where  p  is  the  density, 

q  =  — — o)  •  V<T  (6-2) 

P 

the  potential  vorticity,  oj  the  full  abolute  vorticity,  and  a  the  potential  density  anomaly. 
The  J- vector  in  this  system  is  defined  by 

J  =  pqu  +  F  X  Va - Boj  (6-3) 

9 

in  which  F  encompasses  all  nonconservative  (i.e.,  viscous)  body  forces,  and  B  is  the 
buoyancy  forcing  term, 

I  Da 

B  = - — . 

9  Dt 

In  a  three-dimensional  model,  the  J-vector  comprises  an  advective  component,  a  fric¬ 
tional  term,  and  a  contribution  due  to  non-adiabatic  density  changes. 

It  is  possible  to  demarcate  sub-regions  within  each  layer  and  require  again  that  fluid 
flowing  across  the  equator  between  isentropic  sheets  transform  its  potential  vorticity 
just  as  in  the  shallow-water  model.  As  a  result  there  arises  an  identical  requirement 
to  flux  vorticity  westward  between  different  portions  of  the  boundary  current. 

Since  it  is  the  last  term  in  Equation  6.3  that  distinguishes  a  three-dimensional  from 
a  two-dimensional  model  in  terms  of  the  vorticity  flux,  I  concentrate  on  its  effect  only. 
In  the  absence  of  observational  evidence  that  the  tropics  are  strongly  non-adiabatic, 
it  is  worth  noting  that  in  the  adiabatic  limit,  this  term  vanishes,  and  the  system 
remains  identical  to  the  shallow-water  system.  If  a  buoyancy  flux  exists,  its  ability 
to  flux  vorticity  westward  is  proportional  to  the  buoyancy  itself  and  the  westward 
component  of  absolute  vorticity.  Even  if  the  buoyancy  is  substantial  the  horizontal 
component  of  vorticity  is  quite  small  in  the  absence  of  a  large  vertical  shear.  Thus  it 
is  unlikely  that  this  term  contributes  signiflcantly  to  the  westward  vorticity  flux  in  a 
full  three-dimensional  fluid. 
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This  line  of  reasoning  is  supported  by  the  numerical  experiments  of  Haine  and 
Marshall  (1996)  in  which  the  both  symmetric  and  baroclinic  instability  act  on  a  dif¬ 
ferentially  surface-cooled  fluid  in  an  y^-plane  channel.  In  that  work,  symmetric  insta¬ 
bility  which  has  a  very  fast  time  scale  is  initially  the  primary  mechanism  for  vorticity 
transfer.  However,  baroclinic  instability  ultimately  dominates  the  lateral  transfer  of 
vorticity.  This  process  is  fundamentally  more  efficient  than  symmetric  instability 
and  effectively  shuts  down  the  latter  process.  Symmetric  instability  is  closely  related 
to  the  inertial  instability  that  may  be  present  in  the  tropical  region  in  a  full  three- 
dimensional  fluid,  and  baroclinic  instability  is  quite  similar  to  barotropic  instability. 
Thus  it  seems  reasonable  to  expect  that  although  the  dynamics  may  be  slightly  al¬ 
tered  by  a  local  buoyancy  forcing  in  a  three-dimensional  fluid,  its  contribution  to  the 
vorticity  flux  will  be  small  compared  to  the  eddies  that  emerge  as  a  result  of  the  shear 
instability  of  the  western  boundary  current.  Furthermore,  numerical  experiments  us¬ 
ing  a  three-dimensional  primitive  equation  model  by  Kawase  (1992)  show  basically  the 
same  behavior  observed  in  the  shallow-water  system.  The  boundary  current  becomes 
unstable  at  a  particular  strength  of  the  forcing,  eddies  develop,  and  propagate  across 
the  equator.  There  is  no  indication  that  the  buoyancy  forcing  term  is  significant. 
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Appendix  A 


Numerical  Approximations 


The  time- dependent  numerical  model  is  the  central  tool  by  which  I  analyze  cross- 
equatorial  flow.  Although  the  spatial  discretization  and  issue  of  boundary  conditions 
are  described  in  Chapter  2,  several  important  features  of  the  numerics  are  not  included 
in  that  discussion.  This  appendix  completes  the  description  of  the  model,  clarifying 
the  time-stepping  scheme  used,  the  method  of  stretching  the  model  grid,  and  the 
numerical  approximatins  of  the  vorticity  and  J-vector  components  that  are  essential 
to  the  analysis. 

A.l  Model  Time- Stepping 

The  above  discretization  uses  center-differencing,  and  it  maintains  second  order  spa¬ 
tial  accuracy.  For  time  stepping,  the  essential  feature  is  stability  and  efficiency  over 
accuracy.  To  satisfy  criteria  such  as  the  CFL  condition  and  a  similar  constraint  due 
to  friction,  the  time  step  must  be  quite  short  relative  to  the  time  scale  of  the  motion. 
The  more  severe  limitation  for  numerical  solutions  of  geophysical  systems  in  general 
is  the  spatial  accuracy. 

Arguably,  the  simplest  explicit  time-stepping  scheme  that  is  stable  is  the  three- 
level  modified  Adams-Bashforth  (AB2)  algorithm.  Each  future  state  is  determined 
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by  specific  knowledge  of  the  preceding  two  states  according  to  an  approximation  to 
an  integral  equation.  A  readable  discussion  of  the  original  scheme  and  comparison  to 
other  time  stepping  algorithms  is  found  in  Mesinger  and  Arakawa  (1976). 

The  most  general  time-stepping  problem  is  written  in  the  continuous  form 


i’t  =  $(V’(<)),  (A.l) 

in  which  the  function  $  comprises  all  of  the  terms  that  contribute  to  the  time  rate 
of  change  of  a  variable,  Now  letting  superscripts  denote  time  steps,  t  =  nAt,  the 
discretization  should  approximate 


^n+l 


(A.2) 


The  original  AB2  scheme  approximates  the  right  hand  side  of  Eq.  A.2  as  the  value 
of  $  at  its  midpoint  multiplied  by  At,  as  in  the  leapfrog  scheme,  but  uses  a  Taylor 
expansion  about  the  time  nAt  to  extrapolate  forward  in  time  and  determine  its  value: 


1  /-(n+ljAt 

—  /  ^dt'  ^ 

Ar  JnAt 


-f- 

2  2 


(A.3a) 

(A.3b) 

(A.3c) 


Standard  analysis  of  this  scheme  using  a  simple  oscillation  equation  reveals  that  two 
modes  exist,  one  real  and  one  computational,  as  expected  from  the  fact  that  this 
routine  is  three-level.  The  unphysical  mode  is  damped,  but  the  real  mode  is  unstable, 
for  all  time-steps.  This  feature  is  not  unlike  the  instability  associated  with  two-level 
forward  differencing  (Euler)  in  which  the  term  on  the  right  hand  side  of  Eq  A.2  is 
simply  approximated  by  Just  as  that  instability  is  relieved  by  using  the  implicit 
formulation. 


1  /-(n+ifAt 

At  JnAt 


(A.4) 


the  instability  of  the  AB2  scheme  is  also  relieved  by  shifting  forward  the  time  at  which 
the  forcing  is  to  be  determined.  Thus  the  modified  AB2  scheme  extrapolates  to  a  time 
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slightly  ahead  of  the  mid-point  between  time  steps.  Specifically,  the  time-stepping 
takes  the  form 


5^  =(^e)$--(i +  (A.5) 

where  e  is  a  small  number  (0.1  in  this  study). 

A. 2  Grid  Stretching 

The  original  model  was  specified  on  a  uniform  grid  directed  toward  general  atmospheric 
problems.  Present  knowledge  of  oceanographic  motions  in  confined  basins  allows  an 
improvement  in  the  efficiency  of  the  code  by  stretching  the  grid  and  concentrating 
resolution  in  regions  in  which  the  motion  is  characterized  on  small  spatial  scales. 

Stretching  a  grid  can  be  as  simple  as  varying  grid  dimensions  and  using  the  same 
discretization.  However,  as  discussed  in  Marti  et  al.  (1992),  changing  only  the  grid 
lengths  loses  the  second-order  accuracy  of  the  center-differencing.  Grid  elements  no 
longer  lie  at  midpoints,  where  averaging  and  differential  operations  take  place. 

Finite  difference  formulations  follow  from  Taylor  expansions.  For  the  general  func¬ 
tion,  /(x),  and  letting  primes  denote  derivatives 

f(x  +  Axi)  =  f  (x)  +  Axif  (x)  + 

(*)■■•  ( A.6a) 

/(X-A12)  =  f  (x)  -  Ax2f  (x)  + 

(*)  -  (^)  ■  ■  ■  ■  ( A.6b) 

The  averaging  operator  is  formed  by  adding  Eqs.  A. 6a— A. 6b,  and  the  differential 
operator  by  subtracting  them.  Rearranging  terms, 

/(®)  =  ^(/(a:  + Axi)-l-/(x- Ax2)  + 

(Axi  -  Axz)  f  (x)  +  O  (Axi,  Ax2)„,„^)  (A.7a) 
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n^)  = 


1 


Axi  +  Ax2 


{f{x  +  Axi) 


o((AxO“,(Ax.f)_). 


(A.7b) 


Clearly,  when  increments  of  x  are  equal,  the  schemes  are  of  second  order  accuracy, 
but  not  when  Ai.  are  unequal.  The  only  grid  for  which  all  elements  lie  at  midpoints 
of  adjacent  ones  is  the  uniform  grid.  To  maintain  the  accuracy  of  the  discretization 
on  a  stretched  grid  requires  an  analytical  mapping  of  the  grid.  Following  Marti  et 
al.  (1992),  the  location  of  each  grid  point  is  determined  by  an  analytic  function  (e.g., 
X  —  a:(i))  such  that,  according  to  the  chain  rule, 

d  dx  d 
di  di  dx 


(A.8) 


To  examine  the  accuracy  of  a  scheme  using  this  mapping,  Taylor  expand  the  func¬ 
tion,  g  =  g{i)  about  the  integral  index  i: 


S(i-l)  =  s(i) -5' (•)  +  ■■■  (A.9b) 


As  before, 

SW  =  5(s(i  +  l)  +  s(i-l))  +  0{j")  (A.lOa) 

j'(i)  =  (5(i  +  l)-s(i-l))  +  0(s"')  (A.lOb) 


So  the  mapped  scheme  is  second  order  accurate  in  5'(i).  The  chain  rule  determines 
that 


cPg  fdx\^cPg  (Pxdg 

dP  ^ 

(Pg  f dx\^  <Pg  cPxdx(Pg  (Pxdg 


(A.lla) 

(A.llb) 
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And  it  is  clear  that  the  scheme  is  second  order  accurate  in  g{x)  as  long  as 


O 


O 


O 


(A.  12a) 
(A.12b) 
(A.12c) 


As  long  as  the  change  in  the  variation  of  the  grid  is  sufficiently  small,  the  scheme 
remains  second  order  accurate.  The  argument  extends  straightforwardly  to  more  than 
one  dimension. 

In  the  present  model,  I  chose,  using  %  and  j  as  the  zonal  and  meridional  indices. 


X  =  ax  tanh  (s^  (i  —  to))  +  7x 


(A.13) 


and 

y  =  ay  (tanh  (sj,  {j  -  jo))  +  tanh  (sy  (j  -  ji)))  +  %.  (A. 14) 

The  constants  a^  and  7^  provide  shift  and  scaling  so  as  to  match  the  two  desired 
boundaries,  Si  represent  the  degree  of  stretching,  and  to  and  jo  represent  shifts  for 
the  hyperbolic  tangent  function  to  allow  a  reduction  of  the  grid  separation  near  the 
eastern  and  northern  boundaries  as  well  as  near  the  west  and  equator. 

For  the  runs  presented  below,  using  N  and  M  to  denote  number  of  zonal  and 
meridional  grid  points,  Sx  =  1.85/iV^,  Sy  =  4.3/M,  iO  =  4:Nl5,  jO  =  2M/21.  The 
parameter,  ji  =  M  —  jo,  for  symmetry  about  the  equator.  Figures  A.l  and  A.2 
present  the  grids  used  in  the  standard  runs  of  these  experiments,  with  N  —  255  and 
M  =  511,  along  with  their  derivatives.  Relations  A.12a  -  A. 12c  are  satisfied  by  these 
functions.  Numerical  model  runs  in  Chapter  5  only  use  N  =  127  and  M  =  255,  which 
are  deemed  adequate  to  demonstrate  the  development  of  instability  in  the  vicinity  of 
the  equator. 
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A. 3  Numerical  Approximation  to  Vorticity  Com¬ 

ponents 

Although  the  decomposition  of  the  continuous  system  is  straightforward,  and  the 
vorticity  and  flux  balances  consist  only  of  a  few  terms,  some  errors  result  from  the 
discretization.  The  numerical  algorithm  is  designed  to  conserve  potential  enstrophy, 
not  to  step  forward  the  exact  vorticity  equation.  As  these  errors  are  a  result  of  the 
discretization,  their  magnitude  is  grid-scale  dependent.  In  addition,  they  have  some 
spatial  structure  as  shown  in  a  similar  system  by  Boudra  and  Chassignet  (1988). 

The  numerical  approximation  to  the  vorticity  equation  is 

^  =  Ah  +  SySyC) .  (A.15) 

Errors  are  incurred  because,  for  example, 

4  +  CaS.u^.  (A.16) 

Some  algebra  reveals,  noting  the  commutability  of  the  averaging  and  difference  oper¬ 
ators,  that 

XT’"  IP  s.  (r + r)  +  (r* + f) jiv+ 

(A.17) 

Some  terms  in  Eq.  A.17  are  recognizable  analogs  to  those  of  the  continuous  vorticity 
equation.  However,  some  terms  relating  to  the  horizontal  gradients  of  layer  thickness 
have  no  such  analog.  They  are  simply  a  result  of  the  choice  of  potential  enstrophy 
conservation  in  the  discretization.  As  a  result,  in  the  analysis  of  the  vorticity  budget 
and  J-vectors,  these  errors  must  be  considered. 
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Following  the  method  of  Boudra  and  Chassignet  (1988),  I  decompose  the  vorticity 
equation  into  five  components,  corresponding  to  relative  vorticity  advection,  planetary 
vorticity  advection,  stretching,  friction,  and  the  combined  extra  terms  due  to  the 
conservation  properties  of  the  scheme.  Specifically, 

RVA  =  +  r  (A.i8a) 

pvA  =  (A.i8b) 

STR  =  (A.i8c) 

FRIG  =  -AH[8J^(,-\-8y8y(\  (A.18d) 

EXTRA  =  r' } + 

r'  {«‘"s,r"  +  ji,  [sjff'sj;” + 6~^] }  + 

+ (A.lSe) 

These  are  the  forms  of  relative  vorticity  advection,  planetary  vorticity  advecting, 
stretching,  and  vorticity  associated  with  friction  that  the  model  conserves.  The  EX¬ 
TRA  term  is  simply  the  result  of  the  discretization  of  the  model  and  has  no  physical 
counterpart.  In  the  experiments  RVA,  PVA,  STR,  and  FRIG  are  determined,  and  their 
residual  calculated.  For  all  runs  in  this  study,  the  resolution  of  the  model  is  sufficiently 
high  that  the  residual  is  negligible  compared  to  the  components  that  contribute  to  it. 
All  vorticity  balance  components  are  determined  at  ^-points  of  the  C-grid  model. 

A. 4  Numerical  Approximation  to  J-vectors 

As  with  the  components  of  vorticity,  the  J-vectors  themselves  can  also  be  evaluated 
in  terms  of  the  discretization.  The  model  ensures  that  the  numerical  divergence  of  the 
time  averaged  vector  field 

J  =  -  Ah8^C.  -  AH8yC)  (A.19) 
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vanishes.  As  a  result,  I  define 


(A.20a) 

(A.20b) 

j/„c  =  -A„{s,r‘j.r’)- 

(A.20c) 

The  J-vector  components  lie  at  /i-points  of  the  model. 
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